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ABSTRACT 

A new code for astrophysical magnetohydrodynamics (MHD) is described. The code has been designed to be 
easily extensible for use with static and adaptive mesh refinement. It combines higher-order Godunov methods 
with the constrained transport (CT) technique to enforce the divergence-free constraint on the magnetic field. 
Discretization is based on cell-centered volume-averages for mass, momentum, and energy, and face-centered 
area-averages for the magnetic field. Novel features of the algorithm include (1) a consistent framework for 
computing the time- and edge-averaged electric fields used by CT to evolve the magnetic field from the time- 
and area-averaged Godunov fluxes, (2) the extension to MHD of spatial reconstruction schemes that involve a 
dimensionally-split time advance, and (3) the extension to MHD of two different dimensionally-unsplit inte- 
gration methods. Implementation of the algorithm in both C and Fortran95 is detailed, including strategies for 
parallelization using domain decomposition. Results from a test suite which includes problems in one-, two-, and 
three-dimensions for both hydrodynamics and MHD are given, not only to demonstrate the fidelity of the algo- 
rithms, but also to enable comparisons to other methods. The source code is freely available for download on the 
web. 

Subject headings: hydrodynamics, MHD, methods:numerical 



1. INTRODUCTION 

Numerical methods are essential for the study of a very wide 
range of problems in astrophysical fluid dynamics. As such, the 
development of more accurate and more capable algorithms, 
along with a description of their implementation on modern 
m ' parallel computer systems, is important for progress in the field. 
k> , This paper describes a new code for astrophysical magnetohy- 
j_j ■ drodynamics (MHD) called Athena, developed through a col- 
| laborative effort between the authors. 

There are many numerical algorithms available for solving 
the equations of compressible MHD. One of the most success- 
ful is based on operator splitting of the equations, with higher- 
order upwind methods used for the advection terms, centered- 
differencing for the remaining terms, and artificial viscosity for 
shock capturing. This algorithm, as implemented in for exam- 
ple the ZEUS code (Stone & Norman 1992a; b; Clarke 1996; 
Hayes et al. 2006), has been used for many hundreds of appli- 
cations in astrophysics. The key advantage of the method is its 
simplicity, making it easy to extend with more complex physics 
(for example, Stone & Norman 1992c; Turner & Stone 2001; 
De Villiers & Hawley 2003, Hayes & Norman 2003). 

However, in the fifteen years since the development of 
ZEUS, static and adaptive mesh refinement (SMR and AMR 
respectively) have emerged as powerful techniques to resolve a 

1 Present address: 3915 Rayado PI NW, Albuquerque, NM 87114 



large range in length scales with grid-based methods. Berger 
& Colella (1990) have shown that in order to prevent spurious 
reflections, it is important to enforce conservation at internal 
boundaries between fine and coarse meshes. Thus, operator- 
split methods that do not solve the dynamical equations in con- 
servation form such as ZEUS are unsuitable for use with SMR 
or AMR. This has been our primary motivation for the devel- 
opment of Athena. 

The numerical algorithms in Athena are based on 
directionally-unsplit, higher-order Godunov methods, which 
not only are ideal for use with both SMR and AMR, but also 
are superior for shock capturing and evolving the contact and 
rotational discontinuities that are typical of astrophysical flows. 
Athena is neither the first nor the only MHD code based on 
these methods which is designed for use with AMR; others in- 
clude RIEMANN (Balsara 2000), BATS-R-US (Powell et al. 
1999; Gombosi et al. 2004), AMRVAC (Toth 1996; Nool & 
Keppens 2002), Nirvana (Ziegler 2005), RAMSES (Fromang 
et al. 2006), PLUTO (Mignone et al. 2007), and AstroBEAR 
(Cunningham et al. 2007). While the wealth of papers describ- 
ing AMR MHD codes demonstrates the interest in and impor- 
tance of these numerical methods, it also calls into question 
the need for another paper describing yet another code. How- 
ever, it has been our experience that the precise details of the 
algorithm can be important. The numerical methods in Athena 
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differ, sometimes in small ways, and sometimes in substan- 
tial ways, from those in other codes. Our goals in developing 
Athena have been to write an accurate, easy-to-use, adaptable, 
and maintainable code. Our hope is that the comprehensive de- 
scription provided in this paper will be useful to anyone who 
adopts, modifies, or builds upon the code, as well as for others 
developing their own codes. 

The development of Godunov methods for MHD has re- 
quired substantial progress over the past decade. Most of the ef- 
fort has focused on two main areas: the multidimensional inte- 
gration algorithm, and the method by which the divergence-free 
constraint on the magnetic field is enforced. Different options 
have been explored in different combinations, including uncon- 
strained directionally split integrators (Dai & Woodward 1994), 
or directionally split and unsplit integrators that use either a 
Hodge projection to enforce the constraint (Zachary et al. 1994; 
Ryu et al. 1995; Balsara 1998; Crockett et al. 2005), a non- 
conservative formulation that allows propagation and damping 
of errors in the constraint (Powell 1994; Falle et al. 1998; Pow- 
ell et al. 1999; Dedner et al. 2002), or some form of the con- 
strained transport (CT) algorithm of Evans & Hawley (1988) 
to enforce the constraint (Dai & Woodward 1998; Ryu et al. 
1998; Balsara & Spicer 1999; Toth 2000, hereafter T2000; Pen 
et al. 2003; Londrillo & Del Zanna 2004; Ziegler 2004; Fro- 
mang et al. 2006; Mignone et al. 2007; Cunningham et al. 
2007). T2000 provides a systematic comparison of many of 
these techniques using an extensive test suite. 

While the algorithms in Athena build upon this progress, they 
also incorporate several innovations, including (1) the extension 
of two different directionally unsplit integration algorithms to 
MHD, including the corner transport upwind (CTU) method 
of Colella (1990 - hereafter the CTU+CT algorithm), and a 
simpler predictor-corrector method (see the appendix in Falle 
1981) similar to the MUSCL-Hancock scheme described by 
van Leer (2006; Toro 1999 - hereafter referred to as the VL+CT 
algorithm), (2) the method by which the Godunov fluxes are 
used to calculate the electric fields needed by CT, and (3) the ex- 
tension of the dimensionally-split spatial reconstruction scheme 
in the piecewise parabolic method (PPM) of Colella & Wood- 
ward (1984, hereafter CW) to multidimensional MHD. The 
mathematical foundations of these ingredients for integration 
in two dimensions (2D) is presented in detail in Gardiner & 
Stone (2005, hereafter GS05), and for three dimensions (3D) 
in Gardiner & Stone (2008, hereafter GS08). The focus of this 
paper is on the implementation rather than the mathematics of 
the methods. 

The use of two distinct unsplit integration algorithms in 
Athena, namely the CTU+CT and the VL+CT algorithms, al- 
lows us to compare the advantages and disadvantages of both. 
We find the CTU+CT algorithm is generally less diffusive and 
more accurate than VL+CT. Thus, for simplicity sake, the de- 
scription in this paper will be based on the CTU+CT algorithm. 
However, for some applications the VL+CT algorithm has def- 
inite advantages. A complete description of the 3D VL+CT 
algorithm implemented in Athena, including the results of tests 
in comparison to the CTU+CT algorithm, is provided in a short 



companion paper (Stone & Gardiner 2008, hereafter SG08). 

The primary goal of this paper is to provide a comprehensive 
description of Athena that will serve as a reference for others to 
adopt, modify, and extend the code for their own research. As 
with ZEUS, the source code is freely available from the web, 
along with documentation and an extensive set of test problems 
that are useful for any method. The organization of this paper 
is as follows: §2 introduces the equations of motion solved by 
Athena, while §3 describes their finite-volume and finite-area 
discretizations. Sections 4-6 describe in detail the numerical 
algorithms in one, two, and three spatial dimensions respec- 
tively, including details such as the reconstruction algorithm, 
Riemann solvers used to compute upwind fluxes, and the un- 
split CTU+CT integrator used in multidimensions. In §7 the 
implementation of the algorithms in both C and Fortran95 on 
parallel computer systems is discussed. The results of a com- 
prehensive test suite composed of problems in ID, 2D, and 3D 
are given in §8. Finally, we summarize and discuss future ex- 
tensions to the code in §9. 



2. BASIC EQUATIONS 

Athena implements algorithms which solve the equations of 
ideal MHD, which can be written in conservative form as 



| + V.[pv] = 0, 


(1) 


?^ + V- [pvv-BB+P*l =0, 

ot L J 


(2) 


dE 

— + V- [(£ + P*)v-B(B-v)] =0, 


(3) 


<9B 

— -Vx(vxB) = 0, 
ot 
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where P* is a diagonal tensor with components P* 


= P + B 2 /2 



(with P the gas pressure), E is the total energy density 



and B 2 = B • B. The other symbols have their usual meaning. 
These equations are written in units such that the magnetic per- 
meability /i = 1 . 

An equation of state appropriate to an ideal gas, P= (7— l)e 
(where 7 is the ratio of specific heats, and e is the internal en- 
ergy density), has been assumed in writing equation For a 
barotropic equation of state P = P(p) (for example, P = C 2 p, 
where C is the isothermal sound speed), both equations [3] and [5] 
are dropped from the system. Of course, in this case total en- 
ergy is not conserved. The algorithms implemented in Athena 
can solve the equations of motion in four regimes: both hy- 
drodynamics or MHD with either an ideal-gas or barotropic 
equation of state. In each regime the system of equations to 
be solved is different in number and form, however the same 
general numerical techniques apply. Extension of the numeri- 
cal methods to a more complex, e.g. tabular, equation of state 
is possible. 

It is useful to define vectors of the conserved and primitive 
variables, U and W respectively, with components in Cartesian 
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coordinates (for adiabatic MHD) 



3. DISCRETIZATION 
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where M = p\ is the momentum density. The conservation laws 
can now be written in a compact form (in Cartesian coordinates) 



9U + 9F + aG + 9H_ 

dt dx dy dz 



(7) 



where F, G, and H are vectors of fluxes in the x—, y—, and 
z-directions respectively, with components 



PV X 
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Extension to curvilinear coordinates requires adding metric 
scale factors to the definitions of the fluxes, or using a non- 
conservative formulation that treats grid curvature as source 
terms, or a combination of these approaches. 

For hydrodynamics, or for a barotropic equation of state (or 
for both), the appropriate components of the vectors U, W, and 
their fluxes are dropped. While the last three components of 
these vectors represents the induction equation in Cartesian co- 
ordinates, the numerical algorithm actually used to evolve the 
magnetic field is very different in comparison to that used for 
the other components, as described in the next section. 



Athena integrates the equations of motion on a regular, three- 
dimensional Cartesian grid. The continuous spatial coordinates 
(x,y,z) are discretized into (Nx,N y ,N z ) cells within a finite do- 
main of size (L x ,L y ,L z ) in each direction respectively. The 
cell denoted by indices (i,j,k) is centered at position (xi,yj,Zk)- 
For simplicity we describe the algorithm with the assumption 
that the sizes of the grid cells in each direction, 6x = L x /N x , 
5y = L v /Ny, and Sz = L z /N z respectively, are uniform through- 
out the domain; the numerical methods are easily extended to 
non-uniform grids. 

Time is discretized into N non-uniform steps between the ini- 
tial value fo and the final stopping time tf. Following the usual 
convention, we use a superscript to denote the time level, so 
f" +1 -f" = St". Hereafter we drop the superscript on St with the 
understanding that the time step may vary. 

3.1. Mass, Momentum, and Energy: Finite-Volumes 

Discretizations based on the integral, rather than the differ- 
ential, form of equations Q] through [4] have numerous advan- 
tages for flows that contain shocks and discontinuities (LeV- 
eque 2002). Integration of equation|7]over the volume of a grid 
cell, and over a discrete interval of time St gives, after applica- 
tion of the divergence theorem, 
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is a vector of volume-averaged variables, while 



V(x,y,z,t n ) dx dy dz (12) 
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F(xj-i/2,y,z,t)dydzdt 
(13) 

G(x,yj-i/2,z,t) dxdzdt 
(14) 

H(x,y,z k -i/2,t) dxdydt 



(15) 

are vectors of the time- and area-averaged fluxes. We use the 
convention here, and throughout this paper, that half-integer 
subscripts denote the edges of the computational cells, that is 



is the location of the interface between the cells centered 



at Xi-i and x,-. Thus, the fluxes are evaluated at (and are normal 
to) the faces of each grid cell (see figure 1). Note the half- 
integer superscript on the fluxes denote a time average, rather 
than representing the flux evaluated at f' !+1 / 2 . 

As has been pointed out by many previous authors, equations 
ITTI through [15] are exact: to this point no approximation has 
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been made. A numerical algorithm for MHD within the finite- 
volume approach requires accurate and stable approximations 
for the time- and area-averaged fluxes defined by equations [T3l 
through[l5] In principle, one can approximate the fluxes to any 
order of accuracy, although in practice most algorithms are re- 
stricted to second-order. A variety of authors are exploring the 
use of higher than second-order accurate time- and spatial in- 
tegration (Londrillo & Del Zanna 2000), especially in the con- 
text of WENO schemes (Balsara & Shu 2000; McKinney et al. 
2007). Higher-order schemes improve the accuracy primarily 
in smooth flow, not in shocks or discontinuities, and are more 
difficult to combine with AMR. Based on a set of ID hydro- 
dynamic test problems, Greenough & Rider (2003) conclude 
that a second-order Godunov scheme provides more accuracy 
per computational cost than a fifth-order WENO scheme. Al- 
though it is clear that higher-order schemes will have advan- 
tages for some applications, in Athena we shall restrict our- 
selves to second-order accuracy in both space and time. 

3.2. Magnetic Field: Finite-Areas 

The last three components of equationsQ~T]through[T5]are the 
finite-volume form of the induction equation, which could be 
used to integrate the volume-averaged components of the mag- 
netic field. Instead, in Athena we use an integral form of the 
induction equation that is based on area- rather than volume- 
averages. In GS05, we have argued that area-averaging is the 
most natural representation of the integral form of the induction 
equation. This form conserves the magnetic flux through each 
grid cell, and as a consequence automatically preserves the di- 
vergence free constraint on the field (Evans & Hawley 1988). 

Integration of equation [4] over the three orthogonal faces of 
the cell located at (i— 1 /2,j,k), 1 /2,k) and (i,j,k- 1 /2) 
respectively, gives 
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are the area-averaged components of the magnetic field cen- 
tered on each of these faces, and 
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(24) 

are the components of the electric field £ = -v x B (the elec- 
tromotive force, or emf) averaged along the appropriate line 
element. Note this discretization requires a staggered grid, that 
is the area-averaged components of the magnetic field are lo- 
cated at the faces (not the centers) of the cells. Figure 1 shows 
the relative locations of the cell-centered volume-averaged vari- 
ables (Uijjc), the face-centered area-averaged components of 
the magnetic field (ft,.; i 2j,k, B y,i,j , 2,k,B z ,ij,k , .-) the face- 
centered area-averaged fluxes (F i _ 1 / 2j ^,G / j_ 1 /2, /t ,H /j - it _ 1 /2), 
and the edge-centered line-averaged emfs (£ x ,i,j-i/2,k-i/2> etc.).. 

There are many advantages to using a discretization of the 
induction equation based on area- rather than volume-averages 
(GS05). The most important is that the finite-volume repre- 
sentation, i.e. the cell-volume average, of the divergence-free 
constraint constructed using the time-advanced field 
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(25) 

is kept zero by the discrete form of the induction equation, 
equations [16] through [18] provided of course it was zero at t" 
(Evans & Hawley 1988). Equivalently, the CT algorithm con- 
serves the magnetic flux through each grid cell. The most seri- 
ous disadvantage of using CT with face-centered fields is that it 
complicates the implementation of the algorithm, and the inter- 
face to AMR drivers. 

Of course, there are many possible discretizations of the 
divergence-free constraint, and the CT algorithm based on face- 
centered fields described above preserves only one of them 
(equation [25). T2000 has described an extension to CT which 
preserves the constraint formulated using several different dis- 
cretizations of the divergence operator based on cell-centered 
fields. It is difficult to assess, for a given integration algorithm, 
whether preserving one discretization is more important than 
any other. We have argued (GS05; GS08) that the discretization 
based on face-centered fields is more consistent with the finite 
volume approach in that it conserves the magnetic flux within 
each individual grid cell, equivalently it conserves the volume 
integral of the density of magnetic monopoles at the level of 
grid cells. In addition, in GS08 (see also ® we describe a 
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simple test problem based on the advection of a field loop that 
is sensitive to whether the discretization of the divergence-free 
constraint that is preserved is consistent with the numerical al- 
gorithm used to update the induction equation. If not, growth 
of net magnetic flux will be observed. 

In Athena, the primary description of the magnetic field 
is taken to be the face-centered area-averages equations [T9l 
through [2JJ However, cell-centered values for the field are 
needed to construct the fluxes of momentum and energy (equa- 
tions [8] through [TOl. Here, we adopt the second-order accurate 
averages 



_ 1 

B x ,i,j,k - -j{Bx,M/2,iJi + B x ^\/2,j,k)^ 



By.i.j.k - ~(By.i,j+l/2,k+By.i.j-l/2,k)l 
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(26) 



(27) 



(28) 



Operationally, the face-centered fields are evolved using equa- 
tions [16] through [18] and at the end of each integration step the 
cell-centered fields are computed using equations [26] through 
|28] As shown in GS05 (and discussed further in §5.3), the re- 
lationship between the face- and cell-centered components of 
the field given above determines how their fluxes (the time- and 
line-averaged emfs in equations [22] through [24] and last three 
components of the time- and area-averaged fluxes in equations 
[T3lthrough[T5lrespectivelv) are computed from one another. 

4. ONE-DIMENSIONAL INTEGRATION ALGORITHM 

It is useful (and standard) pedagogy to describe the algorithm 
for integration of the equations of motion in ID first, before in- 
troducing methods for multidimensions. However, for MHD, 
this approach can be misleading. In ID the divergence-free 
constraint reduces to the condition that the longitudinal com- 
ponent of the magnetic field be constant, the CT algorithm is 
not needed, and the discrete forms of the induction equation 
for the area- and volume-averaged fields are identical. As a 
consequence, ID algorithms for MHD are a simple extension 
of those for hydrodynamics. Moreover, ID test problems for 
MHD will not reveal errors associated with the development of 
a non-solenoidal field. Any rigorous test suite for MHD must 
be based on multidimensional problems. 

Nonetheless, we begin a description of the algorithms in 
Athena with the ID integrator as it allows us to introduce ba- 
sic components, such as Riemann solvers and methods for spa- 
tial reconstruction, required in multidimensions. We emphasize 
that the integrators for 2D and 3D MHD, described in detail 
in §5 and §6 respectively, are substantially different and more 
complex than the ID integrator introduced here. 

In ID, the equations of adiabatic MHD can be written in 
Cartesian coordinates as 
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where the vectors of conserved variables and their fluxes are 

pv x 

pv 2 x + P+B 2 /2-B 2 
pv x v y -B x By 

q= M z , f= pv x v z -B x B z , (30) 

(E + P*)v x -(B-v)B x 
B y v x -B x v y 
B z v x — B x v z 

Note these are identical to equations[6]and[8]with the sixth com- 
ponent dropped. We introduce the notation that vectors denoted 
by lower-case letters are in one spatial dimension (and there- 
fore contain 7 components for adiabatic MHD, rather than 8 for 
the same vectors written in full 3D). It is important to remem- 
ber that the components of the ID vectors defined in equation 
l30l will change depending on direction. For example, in the 
y-direction for ideal MHD, the order of the three components 
are permuted (so the second to fourth components become M y , 
M z , and M x respectively), and the sixth and seventh components 
become B z and B x respectively. 

The finite-volume discretization of equation [29] proceeds as 
described in §3.1, giving 

«+l n St / m. 

where 

. 1 
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q; = 



6x 



q(x,t") dx 
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(33) 



is a vector of volume-averaged variables, while 

are the time-averaged fluxes at the interface located at 7c^.\n. 

In a Godunov method, the time-averaged fluxes (equation[7 
are computed using a Riemann solver (see Toro 1999 for an in- 
troduction to the subject). Figure 2 illustrates the process (see 
also LeVeque 2002). Starting from the ID volume-averages 
stored at cell-centers q" a spatial reconstruction scheme is used 
to construct the conserved quantities to the left- and right- 
sides of the interface, q/y-1/2 and qs.,-1/2 respectively. For 
the CTU+CT integrator, the reconstruction is performed in the 
primitive variables, and includes a time-advance using charac- 
teristic variables, with q£,,_i/2 an d <L?,/-i/2 computed from the 
resulting interpolants (this step will be described in detail in 
34.2b . Due to the slope-limiters used to keep the interpolants 
non-oscillatory, the left- and right-states qL,;-i/2 an d Qs.i-i/2 
will not be equal, except in smooth flow. Thus, they define a 
Riemann problem, the solution to which is the time evolution 
of the various waves, and the intermediate states that connect 
them, that propagate away from the interface. The solution to 
the Riemann problem, evaluated at the location of the interface, 
can be used to construct the time-averaged flux (details of the 
calculation of fluxes using Riemann solvers is given in 34.31 ). 



dt dx 



(29) 



4.1. Steps in the ID Algorithm 

The ID algorithm outlined above can be summarized by the 
following steps: 
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Step 1. From q", the volume averages at time level n, com- 
pute the left- and right-states qu-1/2 an d q«,i-i/2 at every in- 
terface using one of the spatial reconstruction algorithms de- 
scribed below in 94.21 

Step 2. Compute the time-averaged fluxes at every interface 

Ci/2 = •7 r (qz,i-i/2,qR : ;-i/2,5 t ,M/2) using one of the Riemann 
solvers described in 94.31 Note the face-centered longitudi- 
nal component of the magnetic field is passed to the Riemann 
solver as a parameter. 

Step 3. Update the cell-centered conserved variables and the 
transverse components of the magnetic field using the finite- 
volume difference equation in ID, equation[3Tl 

Step 4. Increment the time: t n+l =t" + 8t. Compute a new 
timestep that satisfies an estimate of the CFL stability condi- 
tion based on wavespeeds at cell centers 

St = C 6x/max(\v$f \ + C%)) (34) 

where C < 1 is the CFL number, Cp.) is the fast magnetosonic 
speed in the x-direction, evaluated using the updated quanti- 
ties, and the maximum is taken over all grid cells. Note this 
is only an estimate of the CFL stability condition, since the 
wavespeeds used in the Riemann solver can be different from 
those computed from the cell-centered values. 

Step 5. Repeat steps 1-4 until the stopping criterion is 
reached, i.e.. t" +l > tf 

The entire ID integration algorithm is summarized by the 
flow chart shown in figure 3. 

4.2. MHD Interface States 

The first step in the ID algorithm is to compute the left- and 
right-states qu-i/2 an d Qr.;-i/2 tnat define the Riemann prob- 
lem at the interface located at Jt,_i/2. (Note that in our notation 
the left-state qu-i/2 is a ctually on the right side of the cell cen- 
ter at Xi-i, while the right-state q«,,_i/2 is on the left side of the 
cell center at x,-, see figure 2). The reconstruction is inherently 
ID, and therefore is based on the vector of conserved variables 
in ID (equation [30b . This vector contains only the transverse 
components of the field: in ID these are cell-centered quan- 
tities. For reconstruction in multidimensions, the cell-centered 
averages of the face-centered transverse components of the field 
(for example, equations [27] and [28] for reconstruction in the 
x-direction) would be used. When the longitudinal component 
of the field is needed, the area-averaged value stored at the ap- 
propriate interface is adopted. The fact that the longitudinal 
component of the field does not need to be reconstructed from 
cell-centered values is a further advantage of the CT algorithm 
based on staggered (face-centered) fields; it avoids the prob- 
lem of the longitudinal component being discontinuous at the 
interface due to slope-limited reconstruction from cell centers. 

When the CTU+CT unsplit integrator is used in Athena, the 
second- and third-order reconstruction algorithms described be- 
low include both spatial interpolation with slope-limiting in the 
characteristic variables, and a characteristic evolution of the lin- 
earized system in the primitive variables. We have found these 
steps help to make the reconstruction less oscillatory. However, 
they also require an eigenvalue decomposition of the linearized 



equations of motion in the primitive variables. Appendix A cat- 
alogs the eigenvalues and left- and right-eigenvectors for adia- 
batic and isothermal hydrodynamics and MHD in the primitive 
variables needed for this approach. For more complex physics 
(e.g., relativistic MHD) this eigenvalue decomposition may be 
difficult. One advantage of the VL+CT integrator described in 
SG08 is that it does not require a characteristic evolution in 
the reconstruction step. This avoids the need for an eigenvalue 
decomposition in the primitive variables, and therefore this in- 
tegrator may be a better choice for more complex physics. The 
interface state algorithm used in the VL+CT algorithm is de- 
scribed more fully in SG08. 

4.2.1. Piecewise constant (first-order) reconstruction 

The simplest possible reconstruction algorithm is to assume 
the primitive variables are piecewise constant within each cell 
(implying the conserved variables are also piecewise constant), 
leading to the first-order method 

qt,i-i/2 = (35) 

q«,/-i /2 = q< 

First-order reconstruction is far too diffusive for applications, 
however it is useful for testing, or in those circumstances when 
extra diffusion is in fact desired. 

4.2.2. Piecewise linear (second-order) reconstruction 

A better approximation is to assume the primitive variables 
vary linearly within each cell (meaning that the profile of the 
conserved variables within a cell may be steeper than linear). 
This approximation leads to the second-order reconstruction al- 
gorithm used with the CTU+CT unsplit integrator that is given 
by the following steps: 

Step 1. Compute the eigenvalues and eigenvectors of the lin- 
earized equations in the primitive variables using w,, the cell- 
centered primitive variables in ID (which differs from W, de- 
fined in equation[6]only in that it lacks the longitudinal compo- 
nent of the magnetic field). Explicit expressions for these are 
given in Appendix A. 

Step 2. Compute the left-, right-, and centered-differences of 
the primitive variables w, 

Sw u = w/-w/-i, 

SwR,i =W(n-Wi, (36) 
Sw c .i = (w/+i-w,_i)/2 

(Note that in these equations the subscripts L, R, and C refer to 
locations relative to the cell-center at x,.) 

Step 3. Project the left, right, and centered differences onto 
the characteristic variables 

5a L ,i = L(Wi) • 5w L ,i, 

Sa RJ = L(w,) • Sw R:i , (37) 
8a C j = L(Wi) • Syv c ,i 

where L(w,) is a matrix whose rows are the appropriate left- 
eigenvectors computed in Step 1 . 
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Step 4. Apply monotonicity constraints to the differences 
in the characteristic variables, so that the characteristic recon- 
struction is total variation diminishing (TVD), e.g. see LeVeque 
(2002). 

Saf = SIGN( < 5ac,)min(2|<5a L , 1 |,2|5a s ,,|, |«5a c ,,|) (38) 
Step 5. Project the monotonized difference in the character- 
istic variables back onto the primitive variables 

Sw? = 5af • R(wO (39) 
where R(w,) is a matrix whose columns are the appropriate 
right-eigenvectors computed in Step 1 . 

Step 6. Compute the left- and right-interface values using the 
monotonized difference in the primitive variables 

I-max(Af,0)^- 



WL, i+ l/2 = W/ + 



WR,M/2 = Wr 



1 



2Sx 
St 



<5w" 



(40) 



2-^,0)- 



<5w" 



(41) 



where \f and A? are the largest and smallest eigenvalues com- 
puted in Step 1 respectively, at the appropriate cell center. Note 
these values are at different cell faces, with 'Wu+i/i (Wfi,,-i/2) 
located to the right {left) of the cell center at x, . 

Step 7. Perform the characteristic tracing, that is subtract 
from the integral performed in step 6 that part of each wave 
family that does not reach the interface in 8t/2, using (CW; 
Colella 1990) 

St 

Wl.,+1/2 =W L . i+ i /2 + 



\ " 



Wjf,;-i/2 = WRJ-1/2 + 



v; v ^ (( A f- A f)L a -<K , )R a (42) 

X A=>0 

((A?-A?)L a -<SwT)R a (43) 

X A°<0 

where the sums are taken only over those waves that propagate 
towards the interface (i.e., whose eigenvalue has the appropri- 
ate sign), and L" and R Q are the rows and columns of the left- 
and right-eigenmatrices respectively corresponding to A". 

When using approximate Riemann solvers that average over 
intermediate states (like the HLL family of solvers), it is also 
necessary to include a correction for waves which propagate 
away from the interface in order to make the algorithm higher 
than first-order. This is because either the right-interface state 
(if the wavespeed is positive) or the left-interface state (if the 
wave speed is negative) will not include the half-timestep pre- 
dictor evolution in the reconstruction, and will thus be first- 
order. Since the numerical flux in the HLL solver is given by 
a weighted average of the flux in the left-interface state and the 
right-interface state for such waves, the flux itself will be first- 
rder. Specifically, an additional term Aw L (+ i/2 and Aw s ,_i/2 is 
added to each of equations l42land[43lrespectively, where these 
terms are 

Aw L , +1/2 = ~ E ((A? " Af )l a ■ Swf) R" (44) 

A"<0 

Aw s ,,-,/2 = ~ ]T ((A? - A?)L* • <5wf) R a (45) 



A=>0 



We emphasize these terms are not added when the Roe or exact 
solvers are used. 

Step 8. Finally, convert the left- and right-states in the prim- 
itive to the conserved variables, qu-i/2 an d q«,i-i/2- 



4.2.3. Piecewise parabolic (third-order) reconstruction 

Although the numerical algorithms in Athena are formally 
only second-order accurate, we have found that using third- 
order accurate spatial reconstruction can lower the amplitude 
of the truncation error and increase the accuracy of the solution. 
Thus, we have implemented the PPM interface state algorithm 
of CW in Athena. In §|8] we provide a quantitative comparison 
of both the second-order (PLM) and third-order (PPM) recon- 
struction algorithms for smooth and discontinuous solutions in 
ID, 2D and 3D. 

The PPM reconstruction algorithm consists of the following 
steps. 

Steps 1 through 5. These steps are identical to the first five 
steps in the second-order algorithm, see {14.2.21 

Step 6. Use parabolic interpolation to compute values at the 
left- and right-side of each cell center 



y/u = (w, + w ! _ 1 )/2-(<5w;" + ( 5w™ 1 )/6 
w s , ; = (w, +I + w,)/2 - (5<j + (J<)/6 



(46) 



where in the above, the subscript L (R) refers to the left (right) 
side of cell center at x,-. 

Step 7. Apply further monotonicity constraints to ensure 
the values on the left- and right-side of cell center lie between 
neighboring cell-centered values (CW equation 1.10). These 
can be written as a series of conditional statements: 

if (w Ri! -w,)(w, -w t ,,) < 

WL,i = Wj 
Wr,; = W; 

if 6(w Si i-w Lii )(w ! -(w L ,, + WR,,)/2) > (w s .,-w L ,,) 2 

Wt,i = 3w,-2w s ./ 
if 6(w R .,-w L .,)(w, -(w L ,, + w Si! )/2) <-(w Ri ,-w L .,) 2 

w«,,- = 3w,-2wl.; 

These conditions are applied independently to each component 
of w. 

Step 8. Compute the coefficients for the monotonized 
parabolic interpolation function, 

Sw™ = wr.,-w l ,,-, W6,, = 6(w,-(w l ., + w Si! )/2) (47) 

Step 9. Compute the left- and right-interface values using 
monotonized parabolic interpolation (CW equation 1.12) 



WL, i+ l/2 = W Si/ -A n 



2Sx 



St 



Svi" 



1-A n 



W/f,i-l /2 - Wu + A 



,2St\ 



,2St 



^r+(i-A min 3^)w 6 , 



(48) 



(49) 



where A max = max(Af ,0) and A min = min(A?,0) respectively, 
and A^ and Af are the largest and smallest eigenvalues com- 
puted in Step 1 respectively. Note these values are at different 
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cell faces, with w i;+1 /2 (^Rj-1/2) located to the right (left) of 
the cell center at x, . 

Step 10. Perform the characteristic tracing, that is subtract 
from the integral performed in step 9 that part of each wave 
family that does not reach the interface in St/2 (CW; Colella 
1990), using 



w L(i+ i /2 = w i>i+1/2 +^ [l a (A(5yrf-w 6>i )+Bw 6ii )]R a (50) 

A°>0 

w R , +1 /2 = w s , 1+ i/2+5I [L a (C05wr + w 6 , ; ) + Dw 6 ,/)]R Q (51) 



A°<0 



where in the above 



6t_ 

Sx 

~5t_ 

Sx 



(A M A M -A Q A Q ) 



1 2 



(A A - X a X a ) 



where the sums are taken only over those waves that propagate 
towards the interface (i.e., whose eigenvalue has the appropri- 
ate sign), and L" and R a are the rows and columns of the left- 
and right-eigenmatrices respectively corresponding to A". 

Once again, when using the HLL family of solvers, it is nec- 
essary to add a correction for waves which propagate away from 
the interface (as was required in step 7 of the PLM integration). 
These terms are identical to those in equations|44]and|45] which 
are correct to second-order. Again, we emphasize these terms 
are not added when the Roe or exact solvers are used. 

Step 11. Finally, convert the left- and right-states in the prim- 
itive to the conserved variables, qu-1/2 and Q«,!-i/2- 

An important ingredient of the reconstruction algorithm is 
the slope limiters used in steps 4 and 7. It is well-known that 
these limiters clip extrema in the solutions. We have also imple- 
mented the limiters described in Colella & Sekora (2007, here- 
after CS), which are designed to prevent clipping of extrema. 
We find for some tests, the CS limiters significantly improve 
the solution compared to the original PPM limiters used above. 
For the test results shown in §|8]we will always indicate if the 
CS limiters are used. The lesson, however, is that improving the 
convergence rate of the reconstruction algorithm is not always 
the best way to improve the overall accuracy of the solution. 

4.3. Godunov Fluxes 

The second step in the ID algorithm is to compute time- 
averaged fluxes using a Riemann solver. Exact Riemann solvers 
for MHD (e.g. Ryu & Jones 1995) are generally too expensive 
for practical computations with current hardware. Moreover, 
since the full solution to the Riemann problem over all space- 
time is not required, but only the time-integral of the solution 
along the line x = Xf-in (which gives the flux through the inter- 
face), approximate solvers which provide an accurate estimate 
of the flux are all that is needed. In fact, it is not even necessary 
to use the same solver to compute the flux at every interface 
in the grid. Instead, simple solvers can be used in smooth re- 
gions, while more robust (and expensive) solvers are adopted 
only when needed, for example in highly nonlinear flow where 



simple solvers fail (such as strong rarefactions). Since the lat- 
ter generally occupy only a tiny fraction of the total number of 
interfaces over the whole grid, this strategy can be very cost 
effective. 

A wide variety of approximate Riemann solvers for MHD 
are possible, including nonlinear solvers such as the HLL flux 
(Harten et al. 1983), the HLLD flux (Miyoshi & Kusano 2005), 
Toro's FORCE flux (Toro 1999), Roe's linear solver (Roe 1991) 
extended to MHD (Cargo & Gallice 1997), as well as MHD 
solvers based on other approximations (e.g, Dai & Woodward 
1994; 1995; Zachary et al. 1994). A range of solvers is imple- 
mented in Athena, including exact solvers in the simplest cases 
(isothermal hydrodynamics). In the subsections below we de- 
scribe some of the most useful. 

Finally, it is important to emphasize that Godunov methods 
do not require expensive solvers based on complex character- 
istic decompositions. Simple solvers based on the local Lax- 
Friedrichs (LLF) or HLL fluxes that are typically adopted in 
other methods can also be used. Generally, the reason for adopt- 
ing more complex and expensive Riemann solvers is that they 
reduce dissipation, especially in the neighborhood of disconti- 
nuities in the intermediate waves. 

4.3.1. HLL Solvers 

The simplest Riemann solver implemented in Athena uses 
the HLL fluxes as described by Einfeldt et al. (1991), hereafter 
termed the HLLE solver. The HLLE flux at the interface Xp-1/2 
is defined as 



-77 HLLE 
■r 1-1/2 ■ 



b + f L ,i-\/2~b f s ,,-i/2 b + b 

+ ^-^(q«-q«-i) (52) 



b + = max[max(A M , v x ^ R + c R ),0] 
b~ = min [min( A , v x , L - c L ) , 0] 



b + -b~ b + -b 

where f L .,-i/2 = f(qz,,M/2) and f R ,/„i/ 2 = f(qs,i-i/2) are the fluxes 
evaluated using the left- and right-states of the conserved vari- 
ables (using equation[30l), and 

(53) 
(54) 

Here \ M and A are the maximum and minimum eigenvalues 
of Roe's matrix A (see 34.3.21 and Appendix B), v x ,l and v x ,r 
are the velocity component normal to the interface in the left- 
and right-states respectively, and ci and c R are the maximum 
wavespeeds (the fast magnetosonic speed in MHD, or the sound 
speed in hydrodynamics) computed from the left- and right- 
states. The HLLE solver does not require a characteristic de- 
composition of the MHD equations; the eigenvalues of Roe's 
matrix A are given by simple, explicit formulae (see Appendix 
B). Note that if both A M < and v,^ + c R < (or both A > 
and v Xi l — Cl > 0), the HLLE flux will be f^.,-1/2 (or fL,i-i/i)> as 
expected. 

The HLLE solver approximates the solution to the Riemann 
problem using a single constant intermediate state computed 
from a conservative average, bounded using an estimate for the 
maximum and minimum wavespeeds. Thus, for hydrodynam- 
ics it neglects the contact wave, and for MHD it neglects the 
Alfven, slow magnetosonic, and contact waves. For this reason, 
the HLLE is extremely diffusive for these waves (in fact, even 
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if v x = 0, contact discontinuities are diffused with the method). 
Thus, in practice, the HLLE solver is of limited use for appli- 
cations. However, a distinct advantage of the HLLE solver is 
that the intermediate state is positive-definite, that is the pres- 
sure and density in the intermediate state can never be negative. 
Thus, in ID it can be used to construct a positive-definite in- 
tegration algorithm (Einfeldt et al 1991). This is in contrast 
to linearized solvers such as Roe's method, in which the Rie- 
mann solver itself can produce negative densities and pressures 
for one or more of the intermediate states. The HLLE flux is 
therefore an excellent alternative in the rare circumstance that 
a more accurate solver fails. In multidimensions, however, use 
of the HLLE flux at higher than first order does not necessarily 
guarantee the method is positive definite: this depends on the 
details of the multidimensional integrator being used. 

For hydrodynamics, the HLL solver has been extended to in- 
clude the contact wave, resulting in a solution consisting of two 
constant intermediate states bounded by shocks and separated 
by a contact discontinuity. The resulting method is termed the 
HLLC solver. A basic description of the method is given in 
§10.4 of Toro (1999) and will not be repeated here; although 
it is important to note in Athena we choose the wavespeeds 
following the suggestion in Batten et al. (1997). This choice 
has the attractive property that the pressure in the intermediate 
states computed from the Rankine-Hugoniot relations across 
the left and right shocks is the same. We find that for hydro- 
dynamics, this implementation of the HLLC solver produces 
results that are as, if not more, accurate than Roe's method (see 
below), but at much lower computational cost. For ID prob- 
lems, it also is a positive definite method (although again, this is 
not guaranteed in multidimensions). Thus, the HLLC solver is 
highly recommended for adiabatic hydrodynamic simulations 
with Athena. 

Recently, Miyoshi & Kusano (2005) have described an ex- 
tension of the HLL solver to MHD which includes the fast 
magnetosonic, Alfven, and contact waves. The resulting solver 
approximates the solution of the Riemann problem with four 
constant intermediate states. It reduces exactly to the HLLC 
solver when the longitudinal component of the magnetic field 
is zero, and is a positive definite method. The implementa- 
tion of the solver is detailed in Miyoshi & Kusano (2005), and 
will not be repeated here. Tests using Athena indicate that this 
solver, termed HLLD, is typically as accurate as the MHD ex- 
tension of Roe's method, although it is much faster. Thus, the 
HLLD solver is the best choice for many MHD applications 
using Athena. 

4.3.2. Roe's Method 

The HLL fluxes are based on an approximate solution to the 
nonlinear equations of MHD. Instead, Riemann solvers can be 
constructed from exact solutions to an approximate (linearized) 
form of the MHD equations, for example 



dq A/ _,<9q 

57 = A(q) ^- 



(55) 



what makes the system linear). Finding the exact solution to 
linear hyperbolic systems is less difficult because only discon- 
tinuities (no rarefactions) are allowed. 

Of course, the challenge in developing linearized solvers is 
finding the appropriate representation for A(q). Roe (1981) 
proposed one particularly useful linearization, which has sub- 
sequently been extended to adiabatic MHD by Cargo & Gal- 
lice (1997). In this linearization, the Jacobian is evaluated 
using an average state defined in the primitive variables w = 
(p,\,P,B y ,B z ) as follows 



(56) 



v = ( Vft>L + %/7jrv r )/( VpZ+ 

H = (^H L + ^H R )/(^/pT+ Vpr) 

By = (^B y . L + y/pLBy, R )/{y/pE+ y/pl) 

B~z = (VprB z . l + ^/pLB z , r )/(^/pT+ Vpr) 
where H = (E +P*) / p is the enthalpy (used to compute the pres- 
sure), and the subscripts L and R denote the left- and right-states 
of each variable at the interface (computed using one of the re- 
construction schemes described in j34.2t . Explicit forms for the 
matrix A, and its eigenvalues and eigenvectors for isothermal 
and adiabatic hydrodynamics and MHD are given in Appendix 
B. 

Given the eigenvalues A" and left- and right-eigenmatrices 
L(w) and R(w) respectively, where a = 1 ,M denotes the M char- 
acteristics in the solution, the Roe fluxes are simply 



;-l 



/2- 



2 ( fz,,i-l/2+f«,j- 



1 2-rJ2' 



,a I \ a I oa 



(57) 



The matrix A(q) is the Jacobian df/ dq evaluated at some appro- 
priate, constant mean state q (treating this matrix as constant is 



where as before f L ,,-i/ 2 = f(qz.,/-i/2)> ffy-i/2 = f(q«,i-i/ 2 ), and 

a a = L a -S qi . l/2 (58) 

Hi-1/2 = qz.,,-1/2 - qfi,,-i/ 2 (59) 

and the L" and R Q are the rows and columns of the left- and 
right-eigenmatrices corresponding to A a . 

The primary advantage of Roe's method is that it includes all 
of the characteristics in the problem, and therefore is less diffu- 
sive and more accurate than the HLLE solver for intermediate 
waves such as contact discontinuities. Moreover, Roe (1981) 
showed that it gives the flux exactly if the solution to the full 
nonlinear Riemann problem contains only an isolated discon- 
tinuity. However, because it is based on a linearization of the 
MHD equations, for some values of the left- and right-states 
Roe's method will fail (Einfeldt et al. 1991); it will return neg- 
ative densities and/or pressures in one or more of the interme- 
diate states. In Athena, if this occurs we replace the calculation 
of the fluxes at that interface with the HLLE solver (which is a 
positive-definite method) or some other more accurate (e.g. an 
exact) solver. Tests indicate this is only required very rarely. 

5. TWO-DIMENSIONAL INTEGRATION ALGORITHM 

Probably the most popular method for constructing a 2D inte- 
gration algorithm from the ID method described in 50]is based 
on dimensional splitting (Strang 1968). Unfortunately, dimen- 
sional splitting cannot be used for MHD if the equations are to 
be solved in the conservative form. This is because during each 
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one-dimensional update, only the transverse components of the 
magnetic field evolve (e.g., from equation [30l it is clear that B x 
is non-evolutionary during an update in the x-direction). How- 
ever, the divergence-free constraint can only be maintained if 
all three components of the field evolve simultaneously. Thus, 
during the update in the x-direction, B x must evolve. The terms 
that describe this evolution cannot be written in conservative 
form, leading to for example the V • B source term formula- 
tions of Powell (1994) and Powell et al. (1999). However, 
there are significant advantages to maintaining the conservative 
form (T2000), thus in Athena we adopt dimensionally-unsplit 
integrators for MHD, either based on the CTU+CT method (de- 
scribed below), or the VL+CT method (SG08). The use of di- 
rectionally unsplit integrators in multidimensions is one of the 
most important components of the MHD algorithms in Athena. 

Even after adopting an unsplit integration algorithm, combin- 
ing it with the CT method to enforce the divergence-free con- 
straint presents challenges. In particular, the method by which 
the corner-centered, line-averaged emfs are constructed from 
the face-centered, area-averaged fluxes returned by the Rie- 
mann solver is non-trivial. In GS05, we showed that simple 
arithmetic averaging does not work for the unsplit integrators 
adopted here. Instead, we developed several methods for con- 
structing the emfs from the Godunov fluxes, the version actu- 
ally used in Athena is described in 35.31 The resulting method 
reduces exactly to the ID algorithm described in $4] for plane- 
parallel, grid-aligned flow, and preserves the flux normal to the 
plane of the calculation. 

5.1. Steps in the 2D Algorithm 

The 2D CTU+CT integration algorithm is based on the 
method of Colella (1990), and is described in detail in GS05; 
below we provide an overview of the main steps. 

Step 1. Compute and store the left- and right-states at cell 
interfaces in both the x-direction (qL,;-i/2.;,qR,;+i/2.;) and the 
y-direction (qL,u-i/2)<lfi,i,.7+i/2) simultaneously, using any of 
the ID spatial reconstruction algorithms described in 34.21 for 
all the interfaces over the entire grid. Since the ID reconstruc- 
tion algorithms in Athena include a characteristic tracing step, 
when applied in multidimensions the ID reconstruction must 
include V ■ B source terms as described in §3.1 in GS05, and 
briefly in 35.21 Note that the components of (and q«) are 
different on the x— and y-interfaces. 

Step 2. Compute ID fluxes of the conserved variables using 
any one of the Riemann solvers described in 34.3l at interfaces 
in both the x— and y-directions simultaneously 

f/*-l/2j = J 7 (qu-i/2J,qRJ-l/2J,B x j- l /2,j) (60) 
8;*7-l/2 = F(quj-l/2,qRjJ-l/2,By :i j-l/2) (61) 

where the appropriate longitudinal component of the magnetic 
field has been passed to the Riemann solver as a parameter. 

Step 3. Using the algorithm of GS05, described in 35.31 cal- 
culate the emf at cell corners £*,■_! / 2 from the appropriate 
components of the face-centered fluxes returned by the Rie- 
mann solver in step 2, and the z-component of a cell center 
reference electric field £■'" calculated using the initial data at 



Step 4 Evolve the left- and right-states at each interface 
by 5t/2 using transverse flux gradients. For example, the 
mass density, momentum density, energy density, and B z at the 
x-interface located at Xi-in are advanced using 

n+l/2 , St ( „ * \ St 

\t-i/2j = <l£,z-i/2j + (gi-i,;+i/2-gi-W-i/2 J + 2-s A „_(62) 

n+l/2 , St ( „ „ \ St 

Via/ = + 2^ l 8 ''^ 1 / 2 " g «-V2 J + t s,u (63) 

Since the components of ID vectors on the x— and y-interfaces 
differ, care must be taken to associate the components of the 
left- and right-states with the appropriate components of the 
transverse fluxes (for example, the components of qz..;-i/2j with 
the components of g*^ j +l i 2 )- The updates in equations [62] 
and [63] are directionally split (only the transverse flux gradi- 
ent is used) and are based on the conservative form, therefore 
V • B source terms must be added to the momentum density, 
energy, and B z . These are represented by the source term vec- 
tor s A , the last term in both equations. For the left- and right- 
states on the x-interface, the source term vector has compo- 
nents s, : = (0, s M , s E , 0, s B -) where 



„M 



: Bv(/>',.,;.: z,j />',., : 2j) Sx 



s liJ = (BzVz)i,j(B x ,i+l/2,j-B Xl i-i/2j)/Sx 



(64) 



S x'i,j = v z.i.A B x,i+\/2,j-B x ,i-i/2,j)l5x 

Expressions similar to equations [62] and [63] are used to up- 
date the y-interface states located at yj-\n, that is (iu.j-1/2 and 
q s ,j_!/2, for 6t/2 using the flux gradient in the x-direction. 
Source terms analogous to those in equation [64] but propor- 
tional to SBy/Sy, also are necessary (see §4.1.2 in GS05). The 
in-plane components of the magnetic field are evolved using 
CT, 



B' 



,n+l/2 



1 : ~25y~ (^-1/2,7+1/2 _ ^*i-i/2,j-i/2j ( 65 ) 
St 



R n+l/2 _ „ 0t / 

By,iJ-l/2 ~ By,i,j-V2 + ^ ^z.;+l/2j-l/2 ^z,/-l/2,j- 



1/2 



(66) 



using the emfs computed in step 3. 

Step 5. Calculate a cell-centered reference electric field at 
f" +1 / 2 , Si'"* 1 ' 2 , which is needed as a reference state for the CT 
algorithm in step 7. The cell-centered velocities at the half- 
timestep needed to compute £''" +1 ^ 2 come from a conservative 
finite-volume update of the initial mass and momentum density, 
using the fluxes f^Lm; an ^ 8,*>-i/2- ^h e cell-centered compo- 
nents of the magnetic field at the half-timestep come from aver- 
aging the face centered fields at the half-timestep computed by 
equations 1651 and 1661 in step 4 to cell-centers. 

Step 6. Compute new fluxes at cell interfaces using the cor- 
rected left- and right-states from step 4 using one of the Rie- 
mann solvers described in 34.31 giving 



(67) 



time level n, i.e. £ r . 



-(V . B n . —V>. B" ■ ). 



rfH-l/2 _ n+l/2 n+l/2 n n+l/2 s 

Vl/2,; - • / "^Z,,i-l/2,j' < lfi,i-l/2,j'^,i-l/2j- ) 
n+l/2 - / n+l/2 n+l/2 nn+l/2 \ ,^ a -. 
gy-1/2 = T ^%iU/2^RjU^ B y,J-l/2) (68) 

Note the appropriate face-centered fields updated to the half- 
timestep computed in step 4 are passed as parameters to the 
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Riemann solver. If needed, the H-correction is used in this step 
to eliminate the carbuncle instability (see Appendix C). 

Step 7. Apply the algorithm of {35.3l to calculate the CT elec- 
tric fields £**!// 2 j-1/2 usm g tne numerical fluxes from step 6 
and the cell center reference electric field calculated in step 5. 

Step 8. Update the solution from time level n to n+l, us- 
ing the 2D version of the finite-volume difference discretization 
(equation [TT)) for the mass density, momentum density, energy 
density, and B z , and the CT formulae (equations [T6land[T7b for 
the in-plane components of the field B x and B y . 

Step 9. Compute the cell-centered components of the mag- 
netic field from the updated face-centered values using equa- 
tions [26] and [27] 

Step 10. Increment the time: t" +l = t" + 5t. Compute a new 
timestep that satisfies an estimate of the CFL stability condition 
based on wavespeeds at cell centers 



St = C min 



5x 



Sy 



I M+l I ' | v «+l \+C" +1 

I v x,i,j I + fx,i,j I y,i,j I + j , 



(69) 



where C Q < 1 is the CFL number, C^ , and Cfijj are the fast 
magnetosonic speeds in the x- and y-directions respectively, 
evaluated using the updated quantities, and the minimum is 
taken over all grid cells. Note this is only an estimate of the 
CFL stability condition, since the wavespeeds used in the Rie- 
mann solver can be different from those computed from the 
cell-centered values. 

Step 11. Repeat steps 1-10 until the stopping criterion is 
reached, i.e.. t" +l > tf 

The entire 2D integration algorithm is summarized by the 
flow chart shown in figure 4. 

5.2. MHD Interface States in 2D 

In step 1 of the 2D algorithm discussed above, source terms 
must be added to the left- and right-states in the primitive vari- 
ables that arise due to the characteristic tracing step in the re- 
construction algorithms (see ^4.21 ). These terms are necessary 
for a proper accounting of all the evolutionary terms that form 
the characteristic tracing step in multidimensioal MHD (see 
GS05 and GS08 for a complete discussion of the origin of these 
terms). Since the reconstruction is performed in the primitive 
variables, the only terms required are for the transverse compo- 
nents of the magnetic field (in contrast to step 4 in the 2D al- 
gorithm above, where the directional splitting is performed on 
the equations in conservative form, and therefore source terms 
were needed for M, E, and B z ). Thus, for the left-state at the 
x-interface located at X;_i/2, the change to the transverse fields 
due to the source term is 

SBy,u-i/ 2 ,j = j-r,j i., (/>\., i ?., />\., s ?.,) 

while for the left- and right-interface values at the y-interface 
located at yj-i/i, the change to the transverse field due to the 
source term is 



St 



<SB*A«,/-l/2 - ^fy^M- 1 { B y,U-V2-By,i,j-3/2 



Similar expressions are needed for the right-state values at each 
interface (GS05). In both cases the terms are added to the prim- 
itive variables after the reconstruction, and before converting 
back to the conserved variables. 

5.3. Calculating the emfs 

As discussed in $3] the CT update of the magnetic field re- 
quires the line-averaged emfs at cell corners, whereas the Rie- 
mann solver returns area-averaged electric fields at cell faces. 
For example, figure 5 shows the relative positions of the fluxes 
returned by the Riemann solver, and the emfs needed by CT, 
for the 2D grid cell with indices In GS05, it was shown 
that the relationship between the two is determined by the aver- 
aging formulae used to convert between the face-centered area- 
averages of the magnetic field, and the cell-centered volume- 
averages. A variety of different algorithms were explored, 
and the best compromise between accuracy and simplicity was 
found to be 



e. 



-1/2J-1/2 - ^ (^z,i-l/2j + ^,,-l/2j-l+^z,;.>-l/2 + ^,/-lJ-l/2) 




dy / ,_i /2 j-3/4, 
dx 



+ ^11^1 " l IT J C70) 

/-1/4J-1/2 V OX J ,-3/4.^-1/2/ 

where the derivative of £ z on each grid cell face is computed by 
selecting the "upwind" direction according to the contact mode, 
e.g. 



dy 



(dSJdy)^ 
(d£ z /dy)i 



forv t ,,_i /2 >0 
for v,,m/2 < 



i-l/2 



((f),-l + (f),) 0th6rWiSe 



(71) 

(where the subscript j has been suppressed) with an analogous 
expression for the (d£ z /dx). The derivatives of the electric field 
in equation (|7T1 i are computed using the face centered electric 
fields (Godunov fluxes) and a cell center "reference" value £ zi 
e.g. 



d£t 
dy 



J-l/4 



Sy 



(72) 



where the cell center reference electric field ( J is computed 
at the appropriate time level (either t" for step 3 of the 2D algo- 
rithm, or f" +I / 2 for step 7). To help clarify the above, figure 5 
diagrams the relative locations of the Godunov fluxes, corner- 
centered emf, cell-centered reference states, and the derivatives 
of the electric field. Further details are provided in GS05 (and 
GS08 for the 3D case). 

Note for the 3D CTU+CT algorithm, analogous expressions 
to the above are required to convert the x— and y-components 
of the electric field to the appropriate cell corners (see figure 
1). These expressions follow directly from equations[7T|and[72l 
using a cyclic permutation of the (x,y,z) and (i,j,k). 

6. THREE-DIMENSIONAL INTEGRATION ALGORITHM 
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The extension of the dimensionally unsplit CTU integrator 
due to Colella (1990) used in Athena from 2D to 3D is in fact 
quite complex. In particular, for stability with a CFL number 
C < 1 requires 12 Riemann solves per cell per timestep, and 
multiple fractional timesteps are required to correct the left- and 
right-states with transverse flux gradients in a genuinely multi- 
dimensional fashion. This extension of CTU to 3D has been 
described by Saltzman (1994) for hydrodynamics. 

In GS08, we explored the use of the 12-solve CTU+CT al- 
gorithm for MHD, as well as a simpler variant that uses only 
6-solves per timestep, but formally is only stable for CFL num- 
bers C < 0.5. The tests presented in GS08 show that the 6- 
solve algorithm is as accurate as the 12-solve method, and re- 
quires about the same computational cost. However, the 6-solve 
algorithm is dramatically simpler to implement, and therefore 
is the primary 3D integrator used in Athena. 

The 6-solve CTU+CT 3D algorithm is designed in such a 
way that for grid aligned flows it reduces exactly to the 2D 
CTU+CT algorithm described in or the ID algorithm de- 
scribed in $U depending on the symmetry of the problem. Per- 
haps even more importantly, in GS08 we introduced a test prob- 
lem to demonstrate the 3D CTU+CT algorithm preserves a dis- 
crete representation of the divergence-free constraint that pre- 
vents anomalous growth of magnetic flux for problems with 
certain symmetries. The test involves advection of a cylindri- 
cal column of 2D field loops in the x—y plane, with B z = 0, 
and a constant but fully 3D velocity field. In this case the 
z-component of the induction equation reduces to 



i:)B 



dB x ^dBy o 
dt z V dx dy 



Clearly, the second term is proportional to V • B. Thus, if the 
discrete form of the induction equation used to update the field 
components in 3D is able to preserve B Z = Q exactly, then the 
algorithm must preserve the appropriate discrete representation 
of V ■ B = 0. We present results of this field loop test in $8.4l in 
2D, and ES]in 3D. 



6.1. Steps in the 3D Algorithm 

The 6-solve version of the dimensionally unsplit 3D 
CTU+CT algorithm can be described by the following steps 
(see GS08 for details). It may also be useful to compare and 
contrast the steps in the 3D algorithm with those in the 2D 
method (ED. 

Step 1. Compute and store the left- and right-states 
at cell interfaces in the x-direction (qu-i/2,j,k> ( lR,i-i/2.j.k), 
the y-direction ((lL,ij-i/2,ki ( lR,ij-i/2.k)^ an d the z-direction 
( ( lu.j.k-i/2i ( lR,ij,k-i/2) simultaneously, using any of the ID spa- 
tial reconstruction schemes described in $4.21 for all the inter- 
faces over the entire grid. This requires adding V • B source 
terms to the primitive variables, as discussed in GS08 and 36.21 

Step 2. Compute ID fluxes of the conserved variables using 
any one of the Riemann solvers described in $4.3| at interfaces 



in all three dimensions 

V-l/2,;,* = J 7 (qL : i-\/2.jA,<lR,i-\/2J,k,B x j-i/2J,k) (73) 

gy-i/2,* = ^(^mj-iftM^Rjj-i^k^y^j-i^k) (74) 

hy,<t-l/2 = •^'( < lt,i,J,it-l/2 )1fi,i :j,ife-l/2)^z,ij,*-l/2)' (75) 

using the appropriate longitudinal component of the magnetic 
field passed as a parameter to the Riemann solver. 

Step 3. Apply the algorithm of $5.31 to calculate the CT 
electric fields at cell-corners, £*, v _ 1Ait _ 1/2 , £* M/2J ^ 1/2 and 
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, from the appropriate components of the face- 



z,i-l/2,j-l/2,k 

centered fluxes returned by the Riemann solver in step 2, and 
a cell center reference electric field calculated using the initial 
data at time level n, i.e. £-'" k = -(\" jk x B" Jk ). (Note the al- 
gorithms for computing the x— and y-components of the emf 
are a straightforward extension of the algorithm to compute the 
z-component described in $5.31 see GS08.) 

Step 4. Update the face-centered magnetic field by St /2 us- 
ing the CT difference equations [16] through [18] and the emfs 
computed in step 3. 

Step 5. Evolve the left- and right-states at each interface by 
St/ 2 using transverse flux gradients. For example, the hydro- 
dynamic variables (mass, momentum and energy density) are 
advanced using 

n+l/2 _ St / „ t \ 

*iL,i-l/2,j,k ~ HL,i-i/2J,k [Si,j+l/2,k %i,j-l/2,k J 

St f \ St 

~ [ h i.jMi/2 ~ Kik-1/2 ) + y s >.>' i./i (76) 



n+l/2 



-h 

St / 



8i+lj+l/2,* Si+lJ-1/2,* 

St 



k+l/2~K+lJ,k-l/2 ) + y s -t,'j,* (77) 



^R,i-l/2J.k ~ lRJ-l/2J,k 2 c 
6t fh* 

~ 2Sz 

Once again, care must be taken to associate the components of 
the vectors of interface states (e.g. ({u-\/2,j,k) with the appro- 
priate components of the transverse fluxes (e.g., g*^_i/ 2 j an d 
^ijk-1/2)' Moreover, since these updates are directionally split, 
V • B source terms must be added. These are represented by the 
source term vector s x , the last term in both equations. For the 
left- and right-states on the x-interface, the source term vector 
has components s x = (0,s M ,s £ ,0,0) where 



S M ■ -t 

3 x,i,j,k 



dx 



sf jk = -(fiyV^/j^minmod 



where the minmod function is defined as 



dB z 


8B X 


. dz : 


dx 


dB y 


dB x 


dy 


dx 



<J,k 



i.j.k 



minmod (x,y) ■■ 



sign(x)min(|x|, |y|) ifxy>0 







otherwise. 



(78) 



(79) 



The use of the minmod operator to limit the source terms ac- 
cording to the magnitude of the terms in the divergence of B 
is discussed in GS08, it is needed because there are now two 
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terms that arise from transverse gradients, instead of only one as 
in 2D. The transverse components of the magnetic field stored 
at each of the interfaces is evolved using a combination of the 
emfs computed in step 3, and V • B source terms. For example, 
the right-state value of the y— and z-components of the mag- 
netic field at the x-interface at x^.\n are evolved using 

St 



( B y)^-i/2,j,* - ( B y)n,i-i/2j,k- ^ (£*; 

-— (e* 

45z V 



j+l/2,k+l/2 ^x,i,j+l/2,k-l/2 



— f* 

l/2,/t+l/2 c xjJ-l/2,k-l/2 



St. . . JdB z 8B X 



(80) 



^ D <'Rj-l/2J,k ' 



St 



i.j.k 



-a 



: (B z)R,i-l/2 J. k + [^xj,j+l/2,k+l/2~°x,iJ-l/2M+l/2 

St 



4Sy V x '> 



7+1/2,^-1/2 ^x,i,j-l/2M-l/2 



4Sy 

St. . JdB y dB 

-y(^)ij,* mmm Od l 



dx 



(81) 



with similar expressions for the left-state values (but using 
quantities at ; - 1 on the right hand side of the above equations 
as appropriate). The origin of these MHD source terms for the 
transverse components of the magnetic field is discussed fur- 
ther in GS08. The y- and z-interface states are advanced in an 
equivalent manner by cyclic permutation of (x,y,z) and (i,j,k) 
in the above expressions. 

Step 6. Calculate a cell-centered electric field at f" +I / 2 by 
using the fluxes f*_ l/2jr g,*/-i/2,*' and h ;j,*-i/2 t0 com P ut e th e 
cell-centered velocities at the half-timestep using a conserva- 
tive finite volume update for the momentum and density, and 
by averaging the face centered fields at the half-timestep com- 
puted in step 4. This is needed as a reference state for the CT 
algorithm in step 8. 

Step 7. Compute new fluxes at cell interfaces using the cor- 
rected left- and right-states from step 5, and the interface mag- 
netic fields at f" +1 / 2 computed in step 4, using one of the Rie- 
mann solvers described in j 



„«+l/2 

'-l/2j,i : 



n+l/2 



l/2J,jt' ( lR.i-l/2 



«+l/2 



,J,k x,i- 



,n+l/2 



l/2,j,k 



) (82) 



n+l/2 _ T( n+l/2 „" +1 / 2 R n+1/ 2 \ /oo\ 

&i,j-l/2,k ~ ■ ry *iL,iJ-\l2,V*iR,i.j-\l2,V D y,i,j-ll2,k> 



.n+l/2 _ -TV n+l/2 n+l/2 R n+l/2 

.• i m - ■ r \Hmj t k-l/2^R,i,i,k-l/2^ D z,i,i,l 



,) (84) 



ij,jfc-l/2 _ J ^LJJ,k-l/2^R,i,j,k-l/2' 1J z,iJ,k-l/2 J 

using the appropriate longitudinal component of the magnetic 
field passed as a parameter to the Riemann solver. If needed, 
the H-correction is used in this step to eliminate the carbuncle 
instability (see Appendix C). 

Step 8. Apply the algorithm of jj5.3| to calculate the CT elec- 

. • <z i j trn+l/2 c*n+l/2 . s*n+l/2 

trie fields £ x J_ 1/2tk _ 1/2 , £ y!i 4/2J,k-l/2 and S z,i-l/2J-l/2,k usln 8 

the appropriate components of the numerical fluxes from step 7 
and the cell center reference electric field calculated in step 6. 

Step 9. Update the solution from time level n to n + 1 using 
the conservative finite volume update (equation [TT) for the hy- 
drodynamic variables (mass, momentum and energy density) 



and the CT formulae (equations [16] through ITSb to update the 
area-averaged face-centered components of the magnetic field. 

Step 10. Compute the cell-centered components of the mag- 
netic field from the updated face-centered values using equa- 
tions |26]through|28] 

Step 11. Increment the time: t" +l = t" + St. Compute a new 
timestep that satisfies an estimate of the CFL stability condition 
based on wavespeeds at cell centers 

r _ . / Sx Sy Sz 

ot = C mm — 



L,«+l I I r^n+l 

\ V x,i,j,k\~ l ~'-'fx,iJ,k 



v y,i,j,k\ 



\+r n+1 ' \+r n+x 

\ + ^fy,i,j.k \ v zJjM + ^fz,iJ,k, 

(85) 

where C Q < 1 /2 is the CFL number, C" f f ijk , C%[, x . and Cf\ jk 
are the fast magnetosonic speeds in the x— , y—, and z-directions 
respectively, evaluated using the updated quantities, and the 
minimum is taken over all grid cells. Note this is only an esti- 
mate of the CFL stability condition, since the wavespeeds used 
in the Riemann solver can be different from those computed 
from the cell-centered values. 

Step 12. Repeat steps 1-11 until the stopping criterion is 
reached, i.e., t" +l > tf. 

The steps in the 3D integration algorithm are very similar 
to those summarized by the flow chart in figure 4 for the 2D 
algorithm. 

6.2. MHD Interface States in 3D 

As with the 2D integrator, source terms must be added to the 
left- and right-states in the primitive variables calculated using 
the ID spatial reconstruction schemes described in 34.21 Since 
the reconstruction is in the primitive variables, only the trans- 
verse components of the magnetic field require these terms. For 
the right-state at the x-interface located at x,_!/ 2 , the change to 
the transverse fields due to the source terms are 

St { dB dB \ 

(SB y ) R j_ l/2 jM = -— (v y \ M minmod ( -q^~-q^ I (86) 



U,k 



(SB z ) Rti - l/2J M = ~ (v z ),j,*minmod y-g^~-Q^ J ( 87 ) 

Similar expressions are needed for the left-state values, while 
the equations for the left- and right-state values at the y- 
and z-interfaces follow from cyclic permutation of the (x,y,z). 
These terms are added to the primitive variables after recon- 
struction, and before converting back to the conserved vari- 
ables. 

7. IMPLEMENTATION 

The implementation of the numerical algorithms described 
in the previous sections into a functioning computer code can 
be complex, and warrants at least some discussion. 

Athena was developed in C, but many applications scientists 
prefer to work with Fortran. Hence, we have written two dif- 
ferent versions of Athena: the original C code, and another in 
Fortran95. These two versions provide the community with im- 
plementations of the Athena algorithm in the two most popular 
languages used for scientific computing in astrophysics. The 
most important design criteria we have adopted for both ver- 
sions are 
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1. modularity, 

2. documentation, 

3. strict adherence to ANSI standards, 

4. simple control of physics and runtime options 

We briefly discuss each of these below. 

By far the most important design priority is modularity. 
Thus, the Riemann solvers, ID reconstruction algorithms, con- 
version from conserved to primitive variables, boundary con- 
ditions, data output, and the integrators themselves are all bro- 
ken into individual functions, with a common interface specific 
to each class. This makes adding everything from a new Rie- 
mann solver to a new data output format simply a matter of 
writing a new function which conforms to the appropriate in- 
terface. Moreover, all problem-specific code is contained in a 
single file, with functionality provided that makes it easy to add 
new boundary conditions or new source terms in the equations. 

Although writing documentation is never enjoyable, it is crit- 
ical if anyone other than the developer is to use the code. We 
have found this to be true even amongst members of our own 
research groups. The C version of Athena comes with an exten- 
sive User's Guide which describes installing, compiling, and 
running the code, and a Programmer's Guide which explains 
the grid, data structures, and program control and flow. Both 
are included with the source code in the download from the 
web. The Fortran95 version has its own User's Guide. Ample 
comments are also embedded within the source files. 

By adhering to ANSI standards, we ensure Athena can be 
compiled and run on any machine with a C or Fortran95 com- 
piler, as appropriate. To avoid reliance on external libraries, we 
do not use special purpose output formats. The philosophy is 
that data can always be converted into other format by post- 
processing software if needed, or by writing a new user-defined 
output routine. Athena is written to run either as a serial code 
on one CPU or in parallel using domain decomposition through 
MPI calls. The only external libraries needed by Athena are for 
parallelization with MPI (using any version of the MPICH or 
OpenMPI libraries). As algorithms become more complex, the 
use of external libraries for I/O may become unavoidable. For 
example, the HDF5 library has proved to be useful in organiz- 
ing the complex data structures associated with AMR grids. 

The compile and runtime options in the C version of Athena 
are documented in the User's Guide. Physics and algorithm op- 
tions are set at compile time using a configure script generated 
by the autoconf toolkit. In the Fortran95 version, these op- 
tions are determined by selecting which modules to USE. A perl 
build script buildathena is included to simplify the choice 
of problem module, physics, and parallel or serial version. A 
separate user guide is provided for the Fortran code. Both codes 
use a simple block-structured input file with runtime parameter 
values. The Fortran95 version uses NAMELIST and the the C 
version uses a flexible format that emulates NAMELIST func- 
tionality. Although there is nothing special about the specific 
way compiler and run options are set in Athena, the key point 



is that simple and extensible mechanisms to control both are 
provided. 

Two final important aspects of code implementation are the 
single processor performance, and parallelization on distributed 
memory clusters. Aggressive optimization requires mature and 
static algorithms, and often comes at the cost of clarity and 
adaptability in the code. Athena is intended to be a community 
code, and we plan that Athena will continue to be developed 
and extended. Thus, optimization has been limited to the basic 
concepts guided by the rules of data locality and vectorization. 
In the C version, for example, to optimize cache use we define 
all variables within a cell as a data structure, and then create 3D 
arrays of this structure. This ensures values for each variable 
associated with a given cell are contiguous in memory. To pro- 
mote vectorization, as much computational work as is possible 
is done on ID pencils drawn from the grid (for example, the 
spatial reconstruction step). The Fortran95 version is designed 
to take advantage of Fortran array syntax where possible. One 
drawback of dimensionally unsplit algorithms is that the left- 
and right-states and fluxes must be computed and stored for ev- 
ery interface over the entire 3D grid. This requires many 3D 
arrays, which increases the memory footprint of the code and 
reduces cache-performance. However, unsplit algorithms are 
essential for MHD. 

Although Athena requires many more floating point opera- 
tions per cell than algorithms such as ZEUS (as much as ten 
times more), the primary bottleneck on modern processors is 
generally accessing cache and interprocess communication for 
parallel problems. Thus, the performance of Athena in compar- 
ison to ZEUS is not decreased in proportion to the amount of 
work per cell in the two codes. One of the most useful mea- 
sures of performance is the number of cells updated per cpu 
second. This depends on many factors, including the algorithm, 
the size of the grid, and the processor speed. Table 1 lists the 
performance of the C version of Athena on a 2.2 GHz Opteron 
processor, compiled with gcc using an optimization level of 
-03 for various physics and algorithm options and using a 3D 
128 3 grid. For comparison, a 3D version of ZEUS written in 
F77 by one of us (Stone) and run on the same processor gives 
404 x 10 3 cell-updates/sec for adiabatic MHD on a 128 3 grid. 
Thus, while the algorithms in Athena typically require 10 x the 
work of those in ZEUS, the code is only four times slower when 
using the HLLD fluxes. 

Parallelization is achieved in Athena using domain decompo- 
sition with MPI calls to swap data in ghost cells at grid bound- 
aries. The number of ghost cells required depends on the type 
of physics used and the order of the reconstruction. For exam- 
ple, MHD with third-order reconstruction requires four ghost 
cells at every boundary (more are required if the H-correction 
is used, see Appendix C). By sequential exchange of boundary 
conditions in the x—, y—, and z-directions, we avoid the need 
for extra MPI calls to swap values across diagonal domains 
at the corners of the grid. Two factors contribute to making 
Athena very efficient on distributed memory clusters. First, the 
unsplit direct Eulerian update in Athena requires communica- 
tion of ghost zones only once per timestep, greatly reducing the 
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number of MPI calls compared to split methods. Second, the 
ratio of computational work to data communicated is large in 
Athena due to the complexity of the algorithms. Figure 6 plots 
the efficiency of the C version of Athena, defined as the speed 
per processor in a parallel calculation normalized by the speed 
of a single processor calculation, on Red Storm, a Cray XT-3 
at Sandia National Laboratory. Even up to 20,000 processors 
the efficiency of Athena remains above 85%, and is nearly flat 
indicating essentially perfect weak scaling. 

8. TESTS 

Tests are an integral part of the code development process, 
used not only to find bugs in the implementation, but also to 
measure the fidelity of the method in comparison to other tech- 
niques. In this section we present a selection of tests that we 
have found useful in the development of Athena for both hy- 
drodynamics and MHD in ID, 2D and 3D. A more comprehen- 
sive set of tests is published on the web. Many of the problems 
are drawn from test suites of our own codes (Stone et al. 1992) 
or from those published by other authors (Woodward & Colella 
1984, hereafter WC; Ryu & Jones 1995, hereafter RJ; T2000; 
Liska & Wendroff 2003, hereafter LW). Although we begin by 
showing ID tests for hydrodynamics and MHD, our focus will 
be on the multidimensional results that follow, since multidi- 
mensinal tests are so critical for MHD. 

In only a few of the tests do we show the results from more 
than one Riemann solver. In general, we find the most accu- 
rate (and often nearly identical) results are obtained with either 
the Roe and HLLC solvers in hydrodynamics, or the Roe and 
HLLD solver in MHD. Thus, we use these solvers interchange- 
ably. If one solver fails on a particular test, it will be mentioned 
in the discussion. 

8.1. One-Dimensional Hydrodynamics 

Linear wave convergence. One of the simplest, yet most dis- 
criminating tests is to follow the propagation of linear modes 
of each wave family in a periodic domain to measure the am- 
plitude of both diffusion and dispersion errors. Exact eigen- 
functions of sound, contact, and shear waves are initialized in 
a uniform medium with po = 1, Po = 3/5, and 7 = 5/3. The 
wave amplitude A = 10" 6 , and the wavelength is equal to the 
size of the domain L = 1 . For sound waves, the background 
medium is initially at rest. (It is also useful to try a test in which 
v x ,o = ~c s , where c\ = jP/p is the sound speed, so that the right- 
propagating sound waves are standing waves.) For the contact 
and shear waves, the background medium has a constant veloc- 
ity Vxf) = 1 • The solution is then evolved for 1 crossing time, or 
until t/ = l. Figure 7 shows the norm of the Li error vector for 
each wave, defined as 

i 

where q? is the initial solution, as a function of the numerical 
resolution up to 1024 zones, using third-order reconstruction 
and the HLLE, HLLC, or Roe fluxes. The errors for the HLLC 



and Roe fluxes are nearly identical, and converge at second- 
order for each wave family. The errors for the HLLE solver 
are slightly larger, and converge at a slightly lower rate. By 
plotting profiles of the waves, we find the errors are dominated 
primarily by diffusion error; with 16 or more grid points per 
wavelength the plots show almost no dispersion in any of the 
waves. A number of very sensitive tests of the coding can be de- 
signed. Firstly, the Li errors should be identical (to every digit 
of accuracy) for left- and right-propagating waves. Secondly, 
convergence should continue until either the limits of round- 
off error are reached, or nonlinear steeping becomes important 
(when Li ~ A 2 ). We have found that both double precision, and 
very small initial amplitudes, are necessary to see convergence 
out to 1024 cells. This suggests that round-off error can dom- 
inate truncation error in very high resolution simulations with 
higher-order methods such as Athena. 

Sod shocktube. Long a standard test for hydrodynamic codes, 
the Sod shocktube consists of two constant states separated by 
a discontinuity (a Riemann problem). Table 2 lists the values in 
the left- and right-states for this test. Figure 8 shows the results 
for the density, pressure, velocity, and P/p (which is propor- 
tional to the specific internal energy density) at tf = 0.25 when 
run on a grid of 100 cells in the domain -0.5 < x < 0.5 us- 
ing third-order reconstruction, the HLLC Riemann solver, and 
an adiabatic index 7 = 1.4. When configured for ID hydrody- 
namics, Athena reduces to a direct Eulerian PPM code (e.g. §4 
of CW), thus we expect the results should be similar to those 
published by e.g., Greenough & Rider (2003). As is typical of 
a PPM code, Athena resolves the shock front and contact dis- 
continuity with only 2-3 zones. Although we show this test for 
posterity, in our opinion the ID Sod shocktube should no longer 
be considered a discriminating test of algorithms. 

Two- interacting blast waves. Introduced as a test by WC, 
this problem consists of an initially constant density po = 1 in 
a stationary medium in a domain of size L x = 1 with reflect- 
ing boundary conditions, and 7 = 1.4. For x < 0.1, the initial 
pressure is P = 1000, for x > 0.9 P = 100, while P = 0.01 ev- 
erywhere else. The solution is evolved to an arbitrary time of 
tf = 0.038, at which point the shocks and rarefactions gener- 
ated at the two discontinuities in the initial state have interacted 
multiple times in the domain. The test is quite sensitive of the 
ability of the method to capture the interaction of shocks with 
contact discontinuities and rarefactions. Figure 9 shows the 
solution computed with Athena using 400 grid points, third- 
order reconstruction, the CS limiters, and the HLLC Riemann 
solver, with a reference solution computed using 9600 grid 
points shown as a solid line. In addition, the solution can be 
compared to figure 2h of WC. Note that the contact discontinu- 
ity near x = 0.6 is quite smeared out in the Athena solution, this 
seems to be a common property of direct Eulerian methods (see 
figures 18 and 19 in Greenough & Rider 2003), the Lagrange- 
plus-remap version of PPM seems to capture this feature more 
sharply (WC, LW). 

Shu & Osher shocktube. Introduced by Shu & Osher (1989), 
this test measures the ability of a scheme to capture the inter- 
action of shocks with smooth flow. The initial conditions are 
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a strong shock, initially located at x = -0.8, propagating into 
a background medium with a sinusoidally varying density in a 
domain -I <x< \ with adiabatic index 7 = 1.4. Table 2 lists 
the initial conditions for this test. Figure 10 shows the result at 
t = 0.47 computed with both 200 and 800 cells using third-order 
reconstruction, the CS limiters, and the HLLC solver. Compar- 
ison of this plot with, e.g. figure 5 in Balsara & Shu (2000), 
shows the Athena solution is similar to a 3rd-order WENO 
scheme. The use of the CS limiters significantly improves the 
solution in comparison to the original PPM limiters, since with 
only 200 cells many of the extrema in the postshock gas are un- 
resolved, and are clipped with the PPM limiters. We conclude 
that low-order (less than 5th-order) WENO schemes are not 
more accurate than 2nd order Godunov methods like Athena 
for this test. A more comprehensive comparison of Godunov 
and higher-order WENO schemes is provided by Greenough & 
Rider (2003). In particular they conclude for problems involv- 
ing shocks and discontinuities, second-order Godunov schemes 
are more accurate per fixed computational cost. 

Einfeldt strong rarefaction tests. Einfeldt et al. (1991) de- 
scribed several test problems designed to reveal shortcomings 
of various Riemann solvers for hydrodynamics. In particular, 
the Roe solver will always fail on these tests, in the sense that 
it will produce negative densities and pressures in the interme- 
diate states for the initial discontinuity in the first timestep. For 
this reason, when using the Roe solver in Athena we test the 
intermediate states, and if the density or pressure is negative, 
we replace the Roe flux with the HLLE flux for that interface 
only. As an example, figure 1 1 shows the results for the density, 
pressure, velocity, and P/p (which is proportional to the spe- 
cific internal energy) for test 1-2-0-3 in Einfeldt et al. (1991) at 
t = 0.1, computed using 200 grid points, 7 = 1.4, and second- 
order spatial reconstruction (the initial left- and right-states for 
this test are given in Table 2). The profiles of density and pres- 
sure are captured accurately. We find that the HLLE solver is 
only needed for one interface in the first timestep, thereafter the 
Roe solver returns positive states. We have also run the 1-1-2-5 
test in Einfeldt et al. (1991); we find this test is less challenging. 

8.2. One-Dimensional MHD 

Linear wave convergence. As in hydrodynamics, the con- 
vergence of errors in the propagation of linear amplitude MHD 
waves is a sensitive test. For MHD waves, we use a uniform 
medium with p Q = 1, P = 3/5, B = (1 , V2, 1 /2) and 7 = 5/3 in 
a domain of size L = 1 . These choices give well separated wave 
speeds: C/ = 2, Caj = L and C s = 1/2 for the fast, Alfven, and 
slow magnetosonic speeds respectively. Exact eigenfunctions 
for fast and slow magnetosonic, Alfven, and contact waves for 
this background state are given in GS05. These are used to ini- 
tialize each wave family with amplitude A = 10" 6 and exactly 
one wavelength in the domain. Figure 12 shows the norm of the 
Li error vector for each wave family as a function of the numer- 
ical resolution up to 1024 zones, using third-order reconstruc- 
tion and the HLLE, HLLD, or Roe fluxes. The errors using the 
HLLD or Roe fluxes are nearly identical, converge at second- 
order, and are slightly lower than the HLLE fluxes. As before, 



this problem can be used as the basis for a number of very sen- 
sitive tests. For example, standing waves in each family can be 
initialized by setting v A to the appropriate wave speed, the Li 
error should be identical for left- and right-propagating waves, 
and convergence should continue until the limits of round-off 
error or wave-steepening effects are reached. 

Brio & Wu shocktube. An MHD analog to the Sod shock- 
tube was introduced by Brio & Wu (1988), and has now be- 
come a standard test for MHD codes. Table 2 gives the val- 
ues of the primitive variables in the left- and right-states. The 
longitudinal component of the magnetic field is B x = 0.75, and 
is of course constant everywhere. The solution is computed 
with 7 = 2. Figure 13 shows results computed with second- 
order spatial reconstruction and the Roe fluxes, on a grid of 800 
zones at time tf = 0.08. A reference solution, computed using 
10 4 grid points, is shown as a solid line. Once again, shocks 
and contacts are captured in only 2-3 zones. Small oscillations 
are present in the velocity if third-order reconstruction is used, 
indicating our TVD limiters could be improved. Recently, Tor- 
rilhon (2003) has performed a careful study of the convergence 
of finite-volume schemes for MHD Riemann problems similar 
(but not identical) to the Brio & Wu shocktube. We have run 
the regular, nearly coplanar problem defined in §4.2 of that pa- 
per. The left- and right-states for this test are given in Table 
2, in addition B x = 1 . The results, computed using third-order 
reconstruction and the Roe solver, are nearly identical to those 
shown in figure 7 of that paper, although the Athena solution 
with 10 4 grid points is comparable to the solution with twice as 
many points in that paper. At lower resolution (800 grid points) 
the Athena solution shows the compound wave structure which 
appears in dissipative MHD (similar to figure 6 of Torrilhon 
2003). As the numerical resolution is increased, the solution 
converges to the the exact solution for ideal MHD, which does 
not contain this structure. The fact that Athena shows more 
rapid convergence to the exact solution for ideal MHD than the 
central scheme tested in Torrilhon (2003) is indicative of lower 
numerical dissipation. 

RJ shocktube 2a. RJ introduced a large number of MHD 
shocktube problems as tests of a ID algorithm they developed. 
Figure 14 shows the results for the problem shown in their fig- 
ure 2a, which we refer to as the RJ2a test. Table 2 lists the 
left- and right-states for this test, in addition B x = 2. The results 
in figure 14 are computed using third-order reconstruction and 
the Roe fluxes on a grid of 512 cells. This test is of particu- 
lar interest because discontinuities in each MHD wave family 
are produced from the initial conditions, that is both left- and 
right-propagating fast and slow magnetosonic shocks, left- and 
right-propagating rotational discontinuities, and a contact dis- 
continuity. The results in figure 14 show that Athena captures 
each of these discontinuities with 2-4 cells. 

RJ shocktube 4d. A second test introduced by RJ is shown in 
their figure 4d, hereafter we refer to this problem as test RJ4d. 
The left- and right-states are given in table 2, with B x = 0.7. 
The solution at tf = 0.16 is shown in figure 15 computed with 
third-order reconstruction and the HLLD fluxes. The problem 
is interesting because it involves a switch-on slow rarefaction 
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and a slow shock. Although the HLLD solver does not include 
the slow wave explicitly, figure 15 shows these features are all 
captured well in the Athena solution using this solver. 

8.3. Two -Dimensional Hydrodynamics 

Double Mach reflection. Another classic test of hydrody- 
namic algorithms introduced by WC, this problem follows the 
oblique reflection of a Mach 10 shock in air (7 = 1.4). The in- 
teraction of the reflected and incident shocks produces a triple- 
point, and between the resulting contact discontinuity and the 
reflected shock a short jet is formed along the wall. The struc- 
ture of this jet is very sensitive to the numerical diffusion of 
contact waves. This test requires a time-dependent boundary 
condition be applied along the top edge to follow the propa- 
gation of the incident shock; this is easily achieved in Athena 
using function pointers. The problem is initialized following 
the description in WC. Figure 16 shows contour plots of the so- 
lution at t = 0.2 computed with both second- and third-order re- 
construction, and at two different numerical resolutions. The H- 
correction described in Appendix C is used for all the calcula- 
tions to reduce small amplitude noise in the postshock flow. The 
low-resolution (260 x 80) results (first and third panels) show 
small but distinct changes in the jet between the reconstruction 
algorithms. The third-order reconstruction is slightly less dif- 
fusive. Comparison of the results with those in WC (their fig- 
ure 4) demonstrate the differences between the Lagrange-plus- 
remap version of PPM, and the direct Eulerian version imple- 
mented in Athena. The results can also be compared with those 
from ZEUS shown in figures 15 and 16 of Stone & Norman 
(1992a). 

LW implosion test. LW have provided an extensive compar- 
ison of a wide variety of hydrodynamic codes using ID and 
2D problems (including some of the ID problems presented in 
38.lt . We have found the problem discussed in §4.7 in LW, 
hereafter the implosion test, to be one of the most informative. 
It consists of initial states identical to the Sod shocktube prob- 
lem separated by a discontinuity inclined at 45° in a 2D domain 
of size (L x ,L y ) = (0.3,0.3) with reflecting boundary conditions 
everywhere (a more precise description of the initial conditions 
and grid is given in LW). It produces a shock wave which ini- 
tially propagates into the lower left corner, and a rarefaction 
which propagates in the opposite direction. Along the bottom 
and left-side walls, the initial evolution is nearly identical to 
the double Mach reflection test described above. The jets along 
each wall produced in this interaction collide in the lower left 
corner, and produce vortices which propagate outwards along 
the diagonal. In the meantime, a succession of reflected shocks 
interact with the vortices and contact discontinuity, driving the 
Richtmyer-Meshkov instability, and complex shock reflections 
and rarefactions (animations of the evolution, available on the 
Athena web page, are useful for interpreting the evolution). 
Figure 17 shows contours of the density at two times (the same 
two times shown in LW) for a solution computed using third- 
order reconstruction and the HLLC fluxes. The key result of 
the test is the production of the jet along the diagonal. Whether 
this is the correct dynamics was left uncertain in the discussion 



in LW: some codes produce it and others do not. However, we 
have found the jet is reliant on maintaining symmetry in the 
problem. In directionally-split algorithms, perfect symmetry is 
lost, and the collision of the jets in the lower left corner does 
not eject vortices along the diagonal. In dimensionally unsplit 
algorithms such as the CTU method in Athena, the jet is clearly 
formed. We conclude the jet is the correct result, and that it 
is a sensitive test of symmetry. We consider the preservation 
of symmetry a further advantage of the unsplit integrators used 
in Athena, however the primary motivation for their use is the 
preservation of the divergence-free constraint in MHD. 

LW Rayleigh-Taylor instability test. Another test introduced 
by LW in their §4.6 is the nonlinear evolution of a single mode 
of the Rayleigh-Taylor instability. Two fluids, with densities 
two and one respectively, are initialized at rest in a domain of 
size (L x ,L y ) = (1/3, 1) with constant vertical gravitational ac- 
celeration g = 0. 1, and the heavier fluid on top of the light. The 
pressure is computed so that the fluids are in hydrostatic equi- 
librium, with the sound speed equal to one in the light fluid at 
the interface, with 7=1.4. The interface between the two is per- 
turbed with a vertical velocity v y = 0.01 sin(67rx). Running this 
test requires adding gravitational source terms to the equations 
of motion. In Athena, the source terms for a fixed gravitational 
potential are added in such a way as to conserve total energy ex- 
actly, This extension to the algorithms, along with the addition 
of self-gravity in a way that conserves total momentum exactly, 
is described in Gardiner & Stone (in preparation). Without ex- 
plicit viscosity, or surface tension at the interface, there is no 
one correct solution to this problem to which all codes should 
converge. Instead, the resulting structure of the interface be- 
tween the light and heavy fluids is sensitive to the numerical 
diffusion of the method, and to the numerical perturbations in- 
troduced by the grid that seed secondary Kelvin-Helmholtz in- 
stability. Figure 18 shows the results at time tf = 8.5 computed 
with Athena using third-order spatial reconstruction, the HLLC 
fluxes, and a grid of 200 x 400 cells. It can be compared di- 
rectly to the results of other codes shown in figure 4.8 in LW. 
The Athena solution shows more fine-scale structure than many 
other methods, but less than the Lagrange-plus-remap PPM 
codes. This may indicate greater diffusion of contacts in a di- 
rect Eulerian PPM code like Athena, or it may also indicate 
the effect of a contact steepener (which tends to seed more KH 
instability in multidimensions) in the other codes. 

8.4. Two -Dimensional MHD 

Circularly polarized Alfven waves. Circularly polarized 
Alfven waves are an exact nonlinear solution to the equations 
of MHD. T2000 introduced the propagation of these waves as a 
sensitive test of dispersion properties of MHD algorithms. Al- 
though such waves are subject to a parametric instability (Del 
Zanna et al. 2001), for the parameters adopted by T2000 no 
instability should be present. A complete description of this 
test, including the procedure for initializing the solution at an 
oblique angle to the mesh, is presented in GS05. This test 
has proved extremely useful for developing Athena. Figure 19 
shows profiles of the waves after propagating 5 crossing times 
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as a function of resolution, computed using third-order recon- 
struction, the CS limiters, and the Roe fluxes, for both traveling 
and standing waves. Dispersion error is seen to be important 
only at the lowest resolution, diffusion error generally dom- 
inates (this is also true for the linear wave convergence tests 
described in 38. H and 38.2| i. Even with only 8 grid points per 
wavelength, the wave profile is captured well with an ampli- 
tude at least 0.8 of the original. With 16 or more grid points 
per wavelength, the amplitude is better than 0.95 the original in 
both cases. The CS limiter greatly improves the solution at low 
resolution, as it prevents the clipping of extrema in the wave 
profile. Figure 20 shows the norm of the L] error vector as 
a function of resolution for traveling waves, after propagating 
one wavelength, for both second- and third-order reconstruc- 
tion. For comparison, the errors on both a ID and 2D grid are 
shown. In all cases, second-order convergence is evident, with 
the 2D errors larger by a factor of about two. 

Advection of a field loop. This test was introduced and dis- 
cussed extensively in GS05; it consists of the advection of a 
circular field loop by a constant velocity inclined to the grid 
in a periodic 2D domain. For the CT algorithm, solving field 
advection problems is non-trivial. This test demonstrates the 
importance of constructing the line-averaged corner-centered 
emfs used by CT from the area-averaged face-centered elec- 
tric fields returned by the Riemann solver using the technique 
outlined in 35. 3 1 with the CTU integrator. Along with the circu- 
larly polarized Alfven wave test described above, this test has 
been critical to the development of the algorithms. Figure 21 
shows the magnetic field lines and contours of the out-of-plane 
component of the current density J = V x B after advection of 
the loop twice around the domain. The current density is par- 
ticularly sensitive to diffusion or oscillations in the field. The 
figure shows the CTU+CT algorithm in Athena preserves the 
shape of the field loop extremely well. We have also checked 
that if this test is performed with a uniform v z ^ 0, the code 
keeps B z = to round-off error (provided it was zero to begin 
with). As discussed at the beginning of section[6] this confirms 
our formulation of CT preserves the appropriate discretization 
of the divergence-free constraint. 

Orszag-Tang vortex. A 2D MHD test which has now be- 
come a standard is the evolution of the vortex of Orszag & 
Tang (1979). There is some confusion in the literature as to the 
time at which comparisons between solutions are made. The 
results shown here are computed with constant initial densities 
and pressure, po = 25/(367r) and Pq = 5/(1 27r), in a periodic do- 
main of size (L x ,L y ) = (1, 1), with an initial velocity (y x ,v y ) = 
(-sin(27ry), sin(27r:t:)), and a magnetic field computed from the 
vector potential A, = (Bo/47r)cos(47rx)+(.Bo/27r)cos(27ry), with 
Bq = 1/ \/4n. Figure 22 shows contour plots of the density, pres- 
sure, magnetic pressure, and specific kinetic energy density at 
time tf = 1/2 computed on a grid of 192 x 192 cells, which 
can be compared directly to the results in, e.g., T2000 at a time 
of tf = 7r. Of particular note is the symmetry in the solutions. 
Figure 23 shows horizontal slices of the pressure at y = 0.3125 
and y = 0.427 (shown by the horizontal lines in the upper right 
panel of figure 22), with the solution on a 512 2 grid shown as a 



solid line for reference. This test does not seem to be extremely 
discriminating for MHD algorithms. (We consider linear wave 
convergence (see 38.6b . circularly polarized Alfven waves, and 
field loop advection to be more quantitative MHD tests.) The 
most stringent comparison between methods is provided by the 
slices shown in figure 23. Finally, figure 24 plots contours of 
the density, magnetic pressure, specific kinetic energy density, 
and total pressure P* for an isothermal version of the Orszag- 
Tang vortex test. Comparison to results shown previously by 
Balsara (1998, see his figure 8) appear to show significant dif- 
ferences. 

MHD Rotor. The test suite of Stone et al. (1992) contained 
tests based on the propagation of nonlinear amplitude shear 
Alfven waves in ID generated by rotating disks in axisymme- 
try. Since analytic solutions are available for this problem, it 
was possible to provide quantitative measure of the errors in 
ZEUS. (We have confirmed Athena reproduces these tests accu- 
rately, with second-order convergence on the version of the test 
that uses continuous initial conditions.) Following the sugges- 
tion of Brackbill (1986), Balsara & Spicer (1999) introduced 
a 2D version of this test consisting of a rotating disk located 
in the plane of the computation, with an initial magnetic field 
perpendicular to the rotation axis. Strong rotational disconti- 
nuities are generated in the field due to the shear at the sur- 
face of the disk, and shocks and rarefactions are produced by 
the radial expansion of the disk due to unbalanced centrifugal 
forces. We use the initial conditions as described by T2000. 
We present results only for the problem labeled "Rotor Test # 
1", as it involves higher initial velocities and is therefore more 
difficult. No smoothing is used at the surface of the disk. Fig- 
ure 25 shows contours of the density, pressure, Mach number, 
and magnetic pressure at tf = 0.15 on a grid of 400 x 400 cells, 
computed using third-order reconstruction and the Roe fluxes. 
Figure 26 plots slices of the y-component of the magnetic field 
taken along x = 0, and the x-component of the magnetic field 
taken along y = 0. Of note is the near-perfect symmetry main- 
tained in the solutions, with no oscillations. In particular, con- 
tours of the Mach number remain concentric circles in the rar- 
efaction at the center all the way to the origin. Similarly, the 
slices show constant field strength within the central rarefac- 
tion, and sharp discontinuities. 

Magnetic Rayleigh-Taylor instability. To show the effect of 
magnetic fields on the nonlinear evolution of the 2D RT in- 
stability, and to demonstrate the use of AMR with Athena, 
figure 27 shows the results of a single mode RT instability 
computed with 5 levels of refinement. A base grid of 8 x 16 
cells is used, giving an effective resolution on the finest grid 
of 256 x 512. The parameters for this calculation are not iden- 
tical to those used for the LW hydrodynamic RT test shown 
in figure 16. In particular, for the MHD test we use a do- 
main of size (L x ,L y ) = (0.1,0.2) with g = 0.1, an adiabatic index 
7 = 5/3, and densities in the light and heavy fluids of p\ = 1 and 
Ph = 3 respectively. The magnetic field is initially uniform and 
horizontal, with initial amplitude Bo compared to the critical 
value B c = [Lg(ph~ pi)] 1 ^ 2 = 0.14 which suppresses instability 
of Bq/B c = 0.05. The figure shows the distribution of a passive 
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contaminant advected with the flow at a final time tf = 3 in order 
to show mixing, as well as the grid levels used in the AMR cal- 
culation. For reference, the identical calculation but without the 
magnetic field is shown as well. Note the suppression of sec- 
ondary KH instabilities at the interface in the MHD case. An 
extensive discussion of the nonlinear evolution of the magnetic 
RT instability is presented in Stone & Gardiner (2007a; 2007b). 
The use of an AMR grid is very efficient for this problem, since 
the refinement is predominantly required near the interface. 

Blast wave in a strongly magnetized medium. In order to 
demonstrate the propagation of strong MHD shocks in multi- 
dimensions, we show the results of an MHD blast wave prob- 
lem. Many authors have performed similar versions of this test, 
we adopt the initial conditions used in Londrillo & Del Zanna 
(2000). The results are shown at time tf = 0.2 in figure 28 using 
a domain of size (L x ,L y ) = (1,3/2) with a grid of 200 x 300 
cells, third-order reconstruction, and the HLLC (hydro) and 
HLLD (MHD) fluxes. The top row shows contour plots from a 
hydrodynamic version of this test, while the lower row shows 
the MHD results with an initial magnetic field inclined at 45° to 
the grid (B = (Bo/y/2,Bo/y/2) where Bo = 1- By using periodic 
boundary conditions, the flow becomes more complex as the 
outgoing blast wave re-enters the grid on the opposite side, and 
interacts with the contact discontinuity that bounds the evacu- 
ated bubble at the center. Figure 29 shows the result at tf = 1 for 
both the hydrodynamic and MHD problem. Note the CTU in- 
tegrator preserves perfect symmetry (most noticable in the fin- 
gers at the contact discontinuity generated by the Richtmyer- 
Meshkov instability in the unmagnetized problem). Also note 
the magnetic field suppresses the R-M instability (Wheatley et 
al. 2005). Finally, figure 30 plots contours of the MHD blast 
problem using an isothermal equation of state and both Bo = 1 
(top row) and Bo = 10 (bottom row). The plasma (3 = 2P/B 2 = 2 
for Bo = 1, and (3 = 0.02 for Bo = 10 in the external medium 
initially. GS05 shows results for adiabatic MHD with B = 10. 
This problem demonstrates the CTU+CT algorithm is robust 
for low-/3 flows. 

8.5. Three-Dimensional Hydrodynamics 

Noh 's strong shock. As a fully 3D hydrodynamical test, we 
present results from the strong shock test described by Noh 
(1987). This is a very difficult test. A uniform (po = 1), cold 
(P = 0) medium converges in a spherically symmetric radial in- 
flow v r = -1 onto the origin. This generates a very strong (for- 
mally, M = oo) spherical shock wave which propagates away 
from the origin at constant velocity V s = 1/3. Due to the spher- 
ical convergence, the preshock density increases everywhere in 
time according to p(r,t) = po(l + t/r) 2 . However, the density 
immediately upstream of the shock location is always 16, thus 
the postshock gas is uniform with p=64 for 7 = 5/3. A similar 
test is often run in planar (ID) and cylindrical (2D) symmetry, 
however when run with a Cartesian grid the 3D test presented 
here is probably the most difficult. In practice Athena cannot be 
run with pressure identically zero, thus initially we set Po = 10~ 6 
everywhere. The problem is run until tf = 2 in a domain of size 
(L X ,L V ,L Z ) = (1,1,1) computed only in the positive octant with 



200 cells. The inner boundary condition in each dimension is 
reflecting. At the outer boundary the density is evolved accord- 
ing to the analytic solution for the preshock flow, the radial ve- 
locity is held fixed at v r = -1, and the entropy is evolved identi- 
cally to the density, i.e. P(R B , t ) = P (l +r/B B ) 2(1+7) , where R B is 
the spherical radius of the boundary. Figure 3 1 shows contours 
of the density at t = tf computed using second order reconstruc- 
tion, the Roe flux, and the H-correction (see Appendix C). Note 
the contours are smooth and spherical, with virtually no noise 
in the postshock gas. Also shown is a scatter plot of p(r) ver- 
sus r for every eighth grid cell in the computation. The solution 
has the correct density jump and shock speed. The small scatter 
behind the shock demonstrates that with the H-correction, the 
shock remains sharp, smooth, and spherically symmetric. Near 
the origin, the small dip in the density is a signature of wall- 
heating (Noh 1987). These plots can be compared to figure 4.7 
in LW, who ran the same test but in 2D cylindrical symmetry. 
Only a few of the algorithms tested by LW were able to run the 
test at all. The 3D results shown in figure 3 1 are similar to the 
best result in LW (for PPM). Without the H-correction, Athena 
still runs this test but the shock develops strong perturbations 
along the grid directions, similar to but somewhat stronger than 
those evident in the results for the VH-1 code shown in LW. 
Finally, at low resolutions (less than 50 3 ), this test can cause 
Athena to crash when the Roe solver is used, even with the H- 
correction, unless the CFL number is reduced. 

8.6. Three-Dimensional MHD 

Linear wave convergence. We have argued that tests of MHD 
codes must be multidimensional, yet the most quantitative tests 
generally involve plane-wave (ID) solutions. Sensitive multi- 
dimensional tests can be constructed by simply inclining the 
plane wave to the grid at an arbitrary angle. Here, we measure 
the convergence rate of Athena for each MHD wave family in 
3D by initializing a plane waves solution at an oblique angle 
to a grid of size (L x ,L y ,L z ) = (3,3/2,3/2), using the same ini- 
tial conditions as in the ID test described in 38.2l and a grid with 
resolution of2NxNxN cells, with N = 4, 8 , 1 6 , 32 , 64 and 1 28 . 
The angle of the wavevector is chosen so that it does not lie 
along the diagonal of a grid cell, that is there are no symmetries 
in the fluxes in any direction. Details of the initialization of this 
problem in 3D, which requires care to prevent grid noise along 
the wave front, are given in GS08. Figure 32 shows the norm of 
the Li error vector for each wave family using both second- and 
third-order reconstruction computed with the HLLD solver, as 
a function of resolution along L x . For comparison, the errors 
for this same problem in ID are shown as a dashed line. Again, 
we see second order convergence in all wave families. The am- 
plitude of the errors in the fast wave are higher than the ID case 
by about a factor of two, but for all other waves the errors are 
comparable. The fact that the errors in 3D are not significantly 
larger than those in ID reflects the fidelity of the multidimen- 
sional CTU+CT algorithm. 

Circularly polarized Alfven waves. We initialize a ID plane 
wave solution corresponding to the parameter values given by 
T2000 on a grid of size (L x ,L y ,L z ) = (3,3/2,3/2), with the 
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wave front oblique to the grid, following the procedure given 
in GS08. The technique for initializing the wave solution at 
an oblique angle is similar to that used above for linear waves. 
Figure 33 plots profiles of the traveling wave at different res- 
olutions using third-order reconstruction, the CS limiters, and 
the HLLD fluxes. Also shown are the norm of the Li error 
vector computed using both second- and third-order reconstruc- 
tion. These results can be compared directly to the 2D results 
shown in Figure 19. Once again, the solution in 3D compares 
extremely favorably with the 2D solution, for example the Li 
errors are nearly identical to the 2D errors for an adiabatic equa- 
tion of state. 

Advection of afield loop. On a 3D grid, we have found there 
are two challenging versions of this test that can be attempted. 
The first is the 3D analogue of the test described in 38.41 that is 
a cylindrically symmetric field loop with B z = 0, but with a con- 
stant advection velocity along the grid diagonal so that v, y 0. 
As discussed in $6j the numerical algorithm should maintain 
B z = 0, which can only be achieved if the code maintains the bal- 
ance between the two nonzero terms in the z-component of the 
induction equation, that is v z (dB x /dx+dB y /dy) = 0. In turn, for 
constant v z , this requires the code to maintain the divergence- 
free constraint properly. Since the 3D CTU+CT algorithm in 
Athena has been designed to reduce exactly to the 2D version 
for problems with symmetry in z, we obtain the identical re- 
sults for the profile of the field loop in an x—y slice in this test 
as shown in figure 2 1 . Moreover, we confirm that Athena main- 
tains B, = to round-off. A second sensitive test is to incline the 
field loop at an oblique angle to the grid, and advect it with a 
velocity along a perpendicular diagonal (see GS08 for details). 
The resulting current density after advecting the loop twice 
around the grid for both second- and third-order reconstruction 
is shown in Figure 34 for a grid of size (L x ,L y ,L z ) = (1,1,1) 
with 128 3 grid points, and the HLLD fluxes. In this case, the 
component of the field along the axis of the cylinder should re- 
main zero. Although it is not possible to enforce this constraint 
to round-off error (as was the case when the axis of the field 
loop is parallel to a grid direction), nonetheless we find that 
this component is zero to truncation error (see GS08). 

MHD shocktube inclined to the grid. To demonstrate the 
ability of the 3D algorithm to capture shocks and discontinu- 
ities that propagate at an oblique angle to the mesh, we have 
repeated the RJ2a test described in 38.21 on a 3D grid of size 
{L x ,L y ,L z ) = (3/2,1/64,1/64), with the initial discontinuity 
oblique to the grid, using a mesh of 768 x 8 x 8 grid points. 
This gives an effective resolution along the direction of shock 
propagation which is equivalent to the ID test. Initializing the 
discontinuity so as to avoid introducing waves transverse to the 
interface requires care: for more detail see GS08. The results, 
at a time of tj = 1 for the HLLD fluxes and second-order recon- 
struction, are shown in Figure 35. Note that in 3D, each of the 
shocks, contact and rotational discontinuities have been cap- 
tured; there is excellent agreement between the profiles shown 
in figure 35 and the equivalent ID profiles shown in Figure 14. 

Blast wave in a strongly magnetized medium. As a final 3D 
test, we follow the growth of a strong, spherical blast wave in 



a strongly magnetized medium. The initial conditions are iden- 
tical to those given in 38.41 the only difference being that we 
run the problem on a 3D grid of size (L x ,L y ,L z ) = (1, 1.5, 1), 
with 200 x 300 x 200 grid points. Figure 36 shows slices of the 
density and magnetic pressure taken at t = 0.2 computed using 
the HLLD solver and third-order reconstruction. The primary 
difference in the solution compared to 2D is that the size of the 
bubble grows more slowly in 3D, due to the increased adiabatic 
cooling in 3D diverging flow. The contours are all symmet- 
ric and smooth, with no visible asymmetries introduced by the 
grid. 

9. SUMMARY 

We have described Athena, a new code for astrophysical 
MHD. The code implements algorithms based on higher-order 
Godunov methods, with a finite-volume discretization to evolve 
volume-averages of the mass, momentum, and total energy den- 
sity, and a CT algorithm (finite-area) discretization to evolve 
area-averages of the face-centered components of the magnetic 
field. This combination conserves the total mass, momentum, 
energy, and magnetic flux through the grid exactly. Such con- 
servative algorithms are an essential ingredient of AMR meth- 
ods. 

The mathematical foundation of the 2D and 3D algorithms 
in Athena are described more fully in GS05 and GS08. In this 
paper, we have focused on the detailed implementation of the 
methods into a functioning computer code. Step-by-step de- 
scriptions are provided of the multidimensional integrator for 
MHD in 2D and 3D (based on the CTU algorithm of Colella 
1990), the ID reconstruction algorithms (based on an exten- 
sion of the PPM algorithm of CW to multidimensional MHD), 
and a variety of ID Riemann solvers used to compute upwind 
fluxes. We have emphasized the importance of using dimen- 
sionally unsplit integrators for MHD, the advantages of using 
the staggered grid formulation of CT (which requires tech- 
niques for constructing edge-averaged, corner-centered emfs 
from area-averaged face-centered electric fields returned by the 
Riemann solver), and the need to test MHD codes with multi- 
dimensional problems in order to reveal errors associated with 
the divergence-free constraint. 

An extensive series of test problems in ID, 2D, and 3D for 
both hydrodynamics and MHD have been presented. These 
tests, and others published on the web, should be useful to oth- 
ers developing and testing codes for astrophysical MHD. The 
tests show Athena is second-order accurate in space and time 
for smooth solutions in all MHD wave families, even in multi- 
dimensions. We have shown that an advantage of directionally 
unsplit methods is that they preserve symmetries inherent in the 
flow. The 2D CTU+CT method described here reduces identi- 
cally to the ID algorithm for plane-parallel grid-aligned flows. 
Similarly, the 3D CTU+CT method reduces exactly to either 
the 2D or ID methods for plane-parallel, grid-aligned flows, ac- 
cording to the appropriate symmetry. We have exploited such 
symmetries to design a sensitive test of the appropriate stencil 
for maintaining the divergence-free constraint. A planar field 
loop, advected in a fully 3D velocity field, must remain planar. 
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Since the evolution of the component of the field normal to the 
plane of the loop is governed by a term proportional to V ■ B, the 
loop will only remain planar if the divergence-free constraint is 
satisfied exactly on the appropriate stencil. 

In addition to the CTU+CT integrator described in this pa- 
per, an unsplit integrator based on the method described by 
Falle (1981) and similar to the MUSCL-Hancock scheme de- 
scribed by van Leer (2006) has been implemented in Athena. 
The details of this VL+CT method, including tests in 3D and 
comparisons to the CTU+CT method described here, are given 
in SG08. 

The primary motivation for developing Athena has been the 
need to adopt static and adaptive mesh refinement (SMR and 
AMR) to resolve flows over a wide range of length scales in 
various astrophysical applications of interest in our research 
groups (such as magnetized accretion flows, and gravitational 
collapse and fragmentation in dense phases of the ISM) In § 18. 41 
we have shown the results of tests of AMR calculations of 
the Rayleigh-Taylor instability with Athena. Both SMR and 
AMR add considerable complexity to the algorithms, requiring 
special care to conserve mass, momentum, energy, and mag- 
netic flux at fine/coarse grid boundaries. The implementation 
of SMR and AMR with the CTU+CT integrator in Athena will 
be given in a future communication. 

Other extensions to Athena include adding gravitational 
source terms for both a static gravitational potential and self- 
gravity (Gardiner & Stone, in preparation), the shearing box 
(Gardiner & Stone 2006), anisotropic heat conduction (Par- 
rish & Stone 2005; 2007), and transfer of ionizing radiation 
(Krumholz et al. 2007). Many more are either underway (curvi- 
linear coordinates, relativistic MHD, full transport radiation 
MHD), or planned for the future. 



Athena has moved beyond the developmental phase, and is 
now being used for a variety of applications, including stud- 
ies of the MRJ in the shearing box (Gardiner & Stone 2006), 
colliding winds in close binaries (Lemaster et al. 2007), de- 
cay of hydrodynamical turbulence in the shearing box (Shen 
et al. 2006), the magnetic Rayleigh-Taylor instability (Stone 
& Gardiner 2007a; 2007b), shock interactions with magnetized 
clouds (Shin et al. 2007), and the decay of supersonic turbu- 
lence in magnetized molecular clouds (Lemaster & Stone, in 
preparation). 

The Athena code has been made publicly available, and can 
be downloaded from the web, along with extensive documenta- 
tion. Additional test problems beyond those presented here are 
also described on the web. We are confident that Athena will 
become the workhorse for our own applications; it is hoped that 
the description of the algorithms provided in this paper, along 
with the public version of the code provided on the web, will 
be useful to others for solving many problems in astrophysical 
fluid dynamics. 

We have benefited from discussions with many people dur- 
ing the development of Athena; in particular we would like to 
acknowledge Phil Colella, Charles Gammie, Mark Krumholz, 
Nicole Lemaster, Eve Ostriker, and Ian Parrish for their in- 
put. Development of the Athena code was initially supported 
by the NSF ITR program. JS thanks the Royal Society for fi- 
nancial support through the Wolfson Research Merit scheme 
during 2002-2003. Additional support was provided by the 
DOE through DE-FG52-06NA26217. Simulations were per- 
formed on the Teragrid cluster at NCSA, the IBM Blue Gene at 
Princeton University, and on computational facilities supported 
by NSF grant AST-0216105. 
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Table 1 

Performance 1 of Athena in 3D on a 2.2 GHz Opteron processor 



physics 



Roe solver HLLC solver HLLD solver 



isothermal hydro 328 

adiabatic hydro 224 

isothermal MHD 108 

adiabatic MHD 85.9 



340 
242 



124 
97.6 



1 measured in thousands of cell updates per cpu second. 
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Left- and right-states for ID Riemann Problems 
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FIG. 1 . — (left) Centering of volume-averaged conserved variables U and area-averaged components of magnetic field B on the grid, (right) Centering of time- 
and area-averaged components of the fluxes of U, and the time- and line-averaged emfs on the grid. 
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max 



FIG. 2. — {left) An example of piecewise linear reconstruction of conserved variables within each cell to compute the left- and right-states that define a Riemann 
problem at the cell interface. The slopes chosen within each cell are determined by limiters which depend on neighboring cell-center values (not shown) designed 
to prevent the introduction of new extrema. (right) Schematic solution of an MHD Riemann problem in spacetime, consisting of six intermediate states bounded by 
the maximum and minimum wavespeeds. The flux through the interface is the time integral of the solution along the vertical line x = Xi-\n, in this case given by the 
quantities in state qj . In MHD, some characteristics can be degenerate, meaning that the number of intermediate states depends on the problem. 
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FIG. 4. — Flow chart for integration in 2D. The dashed box groups functions that are part of the 2D integrator (described in !5.U . These steps are schematic, with 
many details omitted. The flow chart for the 3D integrator is similar. 
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FIG. 5. — A 2D slice in the x—y plane showing the centering of the fluxes of conserved variables in the x— and y— directions (F and G respectively), and the 
z-component of the emf centered at the cell corner £ z . The CT algorithm used in Athena requires cell-centered reference states for the emf £' to compute the 
gradients (5£/8x) and (5£/5y) which are located between the center of the cell faces and the cell corner. 
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Weak Scaling on Red Storm 
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FIG. 6. — Weak scaling of the efficiency of Athena on a Cray XT-3, using grids with either 32 3 (blue lines), or 64 3 (red lines) cells per processor, and either 
one (SN) or two (VN) processors per node. The quantity \ measures the ratio of the number of cells communicated to the number updated per MPI process. The 
efficiency remains flat well beyond 10 4 processors, indicating excellent weak scaling. 
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FIG. 7. — Convergence in the norm of the Li error vector for sound waves, shear waves associated with each transverse component of velocity, and the entropy 
(contact) wave after propagating a distance of one wavelength in ID. Solutions are computed using third-order spatial reconstruction, and either the Roe fluxes (solid 
line), or HLLE fluxes (dotted line). The errors for solutions computed with the HLLC fluxes are identical to solutions computed with the Roe fluxes. 
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FIG. 8. — Density, pressure, velocity, and specific internal energy (scaled by (7- 1)) for the Sod shocktube test at t = 0.25, computed with 100 grid points, 
third-order spatial reconstruction, and the HLLC fluxes. The solid line is the analytic solution. 
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FIG. 9. — Density at ( = 0.038 in the two-interacting blast wave test, computed with 400 grid points, third-order spatial reconstruction, and the HLLC fluxes. The 
solid line is a reference solution computed with 9600 grid points. 
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FIG. 10. — Density at t = 0.47 in the Shu-Osher Riemann problem, computed with 200 (squares) and 800 (solid line) grid points, third-order spatial reconstruction, 
and the HLLC fluxes. 




FIG. 1 1 . — Density, pressure, velocity, and specific internal energy (scaled by (7- 1)) for the Einfeldt strong rarefaction test at t = 0.1, computed with 200 grid 
points, second-order spatial reconstruction, and the HLLC fluxes. 
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FIG. 12. — Convergence in the norm of the Li error vector for fast, Alfven, slow, and contact waves after propagating a distance of one wavelength in ID. 
Solutions are computed using third-order spatial reconstruction, and either the Roe fluxes (solid line), HLLD fluxes (dashed line), or HLLE fluxes (dotted line). 
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FIG. 13. — Density, pressure, velocity components, transverse component of the magnetic field, and specific internal energy (scaled by (7 — 1)) for the Brio & 
Wu shocktube problem at t = 0.08, computed with 400 grid points, second-order spatial reconstruction, and the Roe fluxes. The solid line is a reference solution 
computed with 10 4 grid points. 
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FIG. 14. — Density, pressure, total energy, all three components of velocity, transverse components and rotation angle <E> = tan l (B z /B y ) of the magnetic field for 
the MHD Riemann problem RJ2a at t = 0.2, computed with 512 grid points, third-order spatial reconstruction, and the Roe fluxes. 
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FIG. 16. — Contours of the density at t = 0.2 for the double Mach reflection test. From top to bottom, the solutions are computed with second order spatial 
reconstruction at low and high resolution, and third order spatial reconstruction at low and high resolution. Here, low resolution uses a grid of 260 X 80 cells, and 
high resolution uses a grid of 520 X 160 cells. All solutions are computed with Roe fluxes and the H-correction. 
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FIG. 17. — Contours of the density at ( = 0.045 (left) and t = 2.5 (right) for the implosion test of Liska & Wendroff. In each case, 31 contours are shown using a 
stepsize of 0.025, starting at a minimum value of 0.125 (at t = 0.045) and 0.35 (at t = 2.5). The solution is computed using third order spatial reconstruction and the 
HLLC fluxes, on a grid of 400 X 400 cells. 
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FIG. 18. — Color image (left) and contours (right) of the density at t = 8.5 in a single mode hydrodynamic Rayleigh-Taylor instability in 2D. Only a single contour 
is shown at p = 1.5 in order to trace the contact discontinuity between the heavy and light fluids. Colors in the image correspond to density values of 0.9 (purple) to 
2.1 (red). The solution is computed using third order spatial reconstruction and the HLLC fluxes on a grid of 200 X 400 cells. 
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FIG. 19. — Profiles of the transverse component of the magnetic field (labelled B2) for both traveling (left) and standing (right) circularly polarized Alfven 
waves, at a time equal to five wave periods, computed on a grid with 2N X N cells, where N = 64 (solid line), 32 (dotted), 16 (dashed line), 8 (dot-dash line), and 4 
(dot-dot-dot-dashed line). Each solution is computed using third order spatial reconstruction and the Roe fluxes. 
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FIG. 20. — Convergence of the norm of the Li error vector for traveling circularly polarized Alfven waves, after propagating a distance equal to one wavelength, 
using an isothermal equation of state. Points marked by squares denote second order spatial reconstruction, triangles denote third order spatial reconstruction. The 
solid lines are solutions computed in ID, the dotted lines are solutions computed in 2D. The dashed line shows the norm of the L| error vector for a 2D solution 
using second order spatial reconstruction computed with an adiabatic equation of state. Also shown is a dashed line with slope -2 for comparison. 
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FIG. 21 . — Magnetic field lines (left) and contours of the z-component of the current density (right) at t = (top row) and alt = 2 after advection of the loop twice 
around the grid (bottom row). The solution is computed using second order spatial reconstruction with the Roe fluxes on a grid of 256 X 128 cells. 
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FIG. 22. — Contours of selected variables at tt = 1/2 in the adiabatic Orszag-Tang vortex test, computed using a grid of 192 X 192 cells, third-order reconstruction, 
and Roe fluxes. Thirty equally spaced contours between the minimum and maximum are used for each plot. The horizontal lines in the panel showing pressure 
correspond to the locations of the slices shown in figure 23. 
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FIG. 23. — Horizontal slices of the pressure at tj = 1/2 in the adiabatic Orszag-Tang vortex test taken at y = 0.3125 (top) and y = 0.427 (bottom). Squares 
correspond to the solution on a 192 X 192 grid, while the solid line is for a 512 2 grid. 
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FIG. 25. — Contours of selected variables at tt = 0.15 in the adiabatic rotor test, computed using a grid of 400 X 400 cells, third-order reconstruction, and Roe 
fluxes. Thirty equally spaced contours between the minimum and maximum are used for each plot. The horizontal and vertical lines in the panel showing magnetic 
pressure correspond to the locations of the slices shown in figure 26. 
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FIG. 26. — Horizontal slice of B y taken at y = (top), and vertical slice of B x taken at x= (bottom) at fy = 0.15 in the rotor test. The solid line is the same data as 
the squares. 
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FIG. 27. — (Left) Grayscale image of the concentration of a passively-advected contaminant at late time in the magnetic Rayleigh-Taylor instability. (Right) Grid 
blocks used to resolve the interface using AMR. The bottom row shows the same quantities, but for a calculation in which the magnetic field strength is zero (i.e., 
hydrodynamics). 
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Fig. 28. — Contours of selected variables at tf =0.2 in the adiabatic blast wave test, computed using a grid of 200 X 300 cells, third-order reconstruction, and either 
HLLC (hydrodynamics, top row) or HLLD (MHD with initial Bo = 1, bottom row) fluxes. Thirty equally spaced contours between the minimum and maximum are 
used for each plot. 
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FIG. 29. — Contours of the density at If = 1 in the hydrodynamic (left) and MHD (right) adiabatic blast test. Fifty equally spaced contours between the minimum 
and maximum are used for each plot. 
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FIG. 30. — Contours of selected variables at tt = 0.2 in the isothermal blast wave test, computed using a grid of 200 X 300 cells, third-order reconstruction, and 
HLLD fluxes. The top row corresponds to an initial Bo = 1. while the bottom row uses an initial So = 10. Thirty equally spaced contours between the minimum and 
maximum are used for each plot. Outgoing waves have already crossed and re-entered the domain by t = 0.2 in the strong field case, thus the contours in the ambient 
medium are due to interaction of these waves rather than oscillations introduced by the algorithm. 




FIG. 32. — Convergence in the norm of the Li error vector for fast, Alfven, slow, and contact waves after propagating a distance of one wavelength at an oblique 
angle across a 3D grid of size 2N X N X N.. Solutions are computed using the HLLD fluxes, and either second-order (solid line) or third-order (dashed line) spatial 
reconstruction. The dotted line shows the errors for second-order spatial reconstruction in ID for reference. 
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FIG. 33. — (Top.) Profiles of the transverse component of the magnetic field for traveling circularly polarized Alfven waves, at a time equal to five wave periods, 
computed on a grid with 2N X N X N cells, where N = 64 (solid line), 32 (dotted), 16 (dashed line), and 8 (dot-dash line). Each solution is computed using third 
order spatial reconstruction and the HLLD fluxes. (Bottom.) Convergence of the norm of the L[ error vector for traveling circularly polarized Alfven waves, after 
propagating a distance equal to one wavelength, for second-order (solid line) and third-order (dashed line) spatial reconstruction. 
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FIG. 34. — Current density in an inclined field loop being advected along the diagonal of a 3D grid at tf = 2 (after twice around the grid). The left panel shows the 
solution for second-order reconstruction, the right for third-order. 




FIG. 35. — Slice through a 3D grid of selected variables for the RJ2a shocktube initialized with the interface oblique to the grid at t = 0.2. This is a fully 3D 
version of the ID test shown in figure 14. 
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FIG. 36. — Contours of selected variables at tf = 0.2 in a 2D slice in the x-y plane at z = (through the center of the grid) in the 3D adiabatic blast wave 
test, computed using a grid of 200 X 300 X 200 cells, third-order reconstruction, and the HLLD fluxes. Thirty equally spaced contours between the minimum and 
maximum are used for each plot. 
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APPENDIX 

EIGENSYSTEMS IN THE PRIMITIVE VARIABLES 



This appendix gives explicit forms for the eigenvalues and eigenvectors of the matrix A resulting from linearizing the dynamical 
equations as W , = A(W)W A , where W is a vector composed of the primitive variables in ID. These eigensystems are needed to 
convert between the primitive and the characteristic variables in the reconstruction algorithms described in jj4.2| 

Adiabatic Hydrodynamics 



(Al) 



For adiabatic hydrodynamics, W = (p, v x , v y ,v z ,P), and the matrix A 
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where a = jP / p {a is the adiabatic sound speed). The five eigenvalues of this matrix in ascending order are 
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Isothermal Hydrodynamics 



For isothermal hydrodynamics, W = (p, v x , v y , v z ), and the matrix A is 



A = 
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v x 








v x 



(A5) 



where C is the isothermal sound speed. The four eigenvalues of this matrix in ascending order are 

A = (v x -C,v x ,v X7 v x + C). (A6) 
The corresponding right-eigenvectors are the columns of the matrix given in equation (A3), with the second column and fifth row 
dropped. The left-eigenvectors are the rows of the matrix 
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Adiabatic Magnetohydrodynamics 
For adiabatic MHD, W = (p,v x ,v y ,v z ,P,b y ,b z ), where b = B/V47T, and the matrix A is 
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where a 2 = jP/p. The seven eigenvalues of this matrix in ascending order are 

A = (v x - C f , v x - Cax, v x - C s ,v x ,v x + Q , v x + C Ax , v x + C f ) 
where the fast- and slow-magnetosonic wave speeds are given by 



Cj, = \ ([a 2 + Cl] ± ^[a 2 + C 2 ] 2 -4a 2 C 2 }j 



(with C/[C S ] given by the +[-] sign). The Alfven speeds are given by 

C 2 = (b 2 x + b] + b 2 z )/p, 

The corresponding right-eigenvectors are the columns of the matrix 
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where S = sign(b x ), and 



c ff = c f a fi 
Qf = C f a f S, 
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In the degenerate case where C A = C Ax = a, so that Cf = C s , then equation (A16) becomes aj = l and a s = 0. The left 
the rows of the matrix 
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(A18) 



(A19) 



are normalization factors for the eigenvectors corresponding to the fast- and slow-magnetosonic waves respectively. 

Isothermal Magnetohydrodynamics 
For isothermal MHD, W = (p, v x , v y ,v z ,by,b z ), where b = B/\/47r, and the matrix A is 
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The six eigenvalues of this matrix in ascending order are 

A = ( Vx - Cf , Vx - C Ax , Vx - Q , Vx + C s , Vx + C Ax , Vx + Cf) , 



(A21) 
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where the fast and slow-magnetosonic wave speeds are given by equation (A10) (with a replaced by the isothermal sound speed C 
here and throughout), and the Alfven speeds are given by equation (Al 1). The corresponding right-eigenvectors are the columns of 
the matrix given in equation (A10), with the fifth row and fourth column dropped. The left-eigenvectors are the rows of the matrix 
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where 
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are normalization factors for the eigenvectors corresponding to the fast- and slow-magnetosonic waves respectively. 
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EIGENSYSTEMS IN THE CONSERVED VARIABLES 

This appendix gives explicit forms for the eigenvalues and eigenvectors of the matrix A resulting from linearizing the dynamical 
equations as U ( = AU^, where U is a vector composed of the conserved variables. These eigensystems are needed to construct the 
fluxes of the conserved quantities using Roe's method (see j H.3.21 i. 

Adiabatic Hydrodynamics 
For adiabatic hydrodynamics, U = (p, pv x , pv y ,pv z ,E), and the matrix A is 



(Bl) 



where the enthalpy H = (E+P)/p, v 2 = v • v, and 7' = (7- 1). The five eigenvalues of this matrix in ascending order are 

A = (v x -a,v x ,v x ,v x ,v x + a), (B2) 

where a 2 = (7- \){H-v 2 /2) = ~fP/p( a is the adiabatic sound speed). The corresponding right-eigenvectors are the columns of the 
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where N a = 1 /(2a 2 ). These are identical to those given by Roe (1981), except the second and third eigenvectors (corresponding to 
the transport of shear motion) have been rescaled to avoid singularities. 

Isothermal Hydrodynamics 
For isothermal hydrodynamics, U = (p, pv x , pv y ,pv z ), and the matrix A is 
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where C is the isothermal sound speed. The four eigenvalues of this matrix in ascending order are 

X = (v x -C,v x ,v x ,v x + C). 
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The corresponding right-eigenvectors are the columns of the matrix given in equation (B3) with the fifth row and fourth column 
dropped, and a replaced by C throughout. The left-eigenvectors are the rows of the matrix 
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where v 2 = v • v, and 
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In these equations 7' = (7- 1), X' = (■y-2)X, Y' = (-y-2)Y, and H = (E + P + b 2 /2)/p. The factors X and Y are introduced by 
the averaging scheme defined by equation (56); thus the matrix A and its eigenvectors depend explicitly on our choice of the Roe 
averaging scheme. The seven eigenvalues of this matrix in ascending order are 

A = ( v x - C f , v x - Cax , v x - C s , v x v x + C s , v x + Cajc , v x + C f ) (B 1 7) 

where the fast and slow-magnetosonic wave speeds are given by 
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The corresponding right-eigenvectors are the columns of the matrix 
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where the C// iJ . s ,g/,.s J A/ i . s ,a/ iS and P ytZ are given by equation (A13) to (A17) (with a replaced by a), V,y jS = viaj^ (i = x,y,z), and 

^51 = a, (H'-v x C f ) + QAvyp; + v z (i* z )+A s blP* 2 / Pl (B22) 
#52 = -(v y P z - v z fi y ) = -R 56 , (B23) 
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R 5 , = a s (H' - v x C s ) - Qf(v y (3* + v z fi*) -A fb* ± P* 2 / p, 

R 54 = v 2 /2+X'/ 1 ' 
R 55 = a s (H' + v x C s ) + Q f (v y (3; + v z (3:)-A f blPl 2 /p, 
R 57 = a f (H' + v x C f ) - QAvyp; + v z /3*)+A s bl/3? / p. 



where H 1 =H-b 2 / p. In these equations 
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The left-eigenvectors are the rows of the matrix 












■ L n -V xf -C ff -V yf + Q s Q* y 


-%f+QsQt 


Olf 


A s Q*-a f b y 
-frS^p/2 


A s Q*-a f b z 




La ^ -p z /2 


Py/2 





PySJp/2 




L31 -V xs -C ss -V ys -QfQ* 


-%-QfQZ 


a s 


-A f Q*-a s b y 


-A f Q*-a s b z 


L = 


L41 2v x 2v y 


2v z 


-i/a 


2 2by 


2b z 




L51 -V xs +C ss -V ys + Q f Q*y 


-%+QsQ* 


a s 


-A f Q*-a s by 


-A fQ* - a s b z 




L 6 i p z /2 


-py/2 





-PzSJp/2 


PyS^p/2 ' 




_ hn -V xf +C ff -V yf -Q S Q; 


-%f-QsQ: 


a f 


A s Q*-a f b y 


A s Q*-a f b z _ 



where a symbol over the quantity q denotes normalization via q = Yq/(2a 2 ) orq = q/(2a 2 ). In addition, 

q; = p;/p?, q;=p;/p?, 

and 

L n = a f (v 2 - H') + C ff (C f + v x ) - Q s {v y Q* y + v z Q\) -A,\b± \ /p, 

L2i=(v y (3 z -v z (3y)/2 = -L 6i 
L 31 = a s (v 2 -H') + C ss (C s + v x ) + Q f (vyQ; + v z Q;)+A f \b ± \/p, 

L 4l = l-v 2 + 2X' 
L51 = a s (v 2 -H') + C ss (C s -v x )-Q f (vyQ; + v z Q:)+A f \b ± \/p, 
L 7l = a f (v 2 - H') + C ff (C f - v x ) + Q s (v y Q; + v z Q*) -A,\b± \ / p, 

Isothermal Magnetohydrodynamics 
For isothermal MHD, U = (p, pv x , pv y , pv z ,b y ,b z ), where b = B/\/47r, and the matrix A is 






1 











■ 


-v 2 x + C 2 +X 


2v x 








b y Y 


b z Y 


~V X Vy 


Vy 


Vx 





-b x 





-v x v z 


v z 










-b x 


(b X Vy-byV X )/p 


by/P 


-b x /p 





Vx 





(b x v z -b z v x )/p 


b z /p 





-bx/p 





Vx . 



(B24) 

(B25) 
(B26) 
(B27) 

(B28) 



(B29) 



(B30) 

(B31) 
(B32) 
(B33) 
(B34) 
(B35) 

(B36) 



(B37) 



where C is the isothermal sound speed, and X and Y are given by equations (B15) and (B16). The six eigenvalues of this matrix in 
ascending order are 

X = (v x -C f ,v x - C Ax , v x - C s ,v x + C s ,v x + C Ax , v x + Cf) (B38) 
where the fast- and slow-magnetosonic wave speeds are given by 

^2 



-f,s ■ 



\([c 2 +c 2 A ] 



± 



c 2 +c 2 ] 2 - 



*c 2 cj x 



(with Cf[C s ] given by the +[-] sign), where C 2 = C 2 +X, and the Alfven speeds are 



C l A =C z Ax + b* ± 2 /p 



C Ax = K/p, 



b* 2 = Y(b z y + bp. 



(B39) 



(B40) 



The corresponding right-eigenvectors are the columns of the matrix given by equation (B21) with the fifth row and fourth column 
dropped, and a replaced by C in the definitions given in equations (A15)-(A16). The left-eigenvectors are the rows of the matrix 



(B41) 



Ln 


~ c ff 


QsQ; 

-Pz/2 


QsQ: 


A S Q*y 

-PzSy/p/2 

-AfQy 


Asq: 


(v y Pz-v z (3 Y )/2 





Py/2 


PyS^p/2 

-AfQ* 




-C 


-QfQy 


-QfQ*z 




c 

ss 


QfQ*y 


Q/Q* 


-AfQ*y 


-A f Q* z 


-(VyP Z -V z Py)/2 





Pz/2 


-Py/2 


-PzSVp/2 


PyS^p/2 


L(,\ 


c ff 




-&Q* 


a s q; 


AsQ* z _ 



61 



where Cff. ss ,Qf. s , and A/ s are given by equations (A13)-(A15) (with a replaced by C), (3 y . z are given by equation (A17), Q* are 
given by equation (B30), and 

Ln=Cff(C f + v x )-Q s (v v Q; + v z Q:)-A s \b* ± \/p, (B42) 

L 31 = CJC s + v x ) + Q f (v y Q; + v z Q;)+A f \bl\/p, (B43) 

L 4 i = C^CC-vJ-Q/Cv^ + v.ep+A/lfoll/p, (B44) 

L 61 =C // (C / -vJ + 4(v y G;+v,e z *)-A>l|/p. (B45) 
In these equations, a symbol over the quantity q denotes normalization via q = q/(C 2 [l + «/]) and q = q/(C 2 [l + a 2 ]). 

THE H-CORRECTION: FIXING THE CARBUNCLE PROBLEM 

For strong, planar shocks in multidimensions propagating along a grid direction, higher-order Godunov methods can be subject to a 
numerical instability (Quirk 1994) that grows into large amplitude perturbations of the shock front at the grid scale. This "carbuncle" 
instability can easily be demonstrated with a simple 2D test: a uniform high Mach number flow in the +x-direction is initialized 
everywhere in the domain, with inflow boundary conditions on the right boundary, and reflecting everywhere else. If zone-to-zone 
perturbations in the density with small amplitude (Sp/p= 10~ 4 ) are added, the reflected shock will develop the carbuncle instability 
as it propagates to the left across the grid. Radiative cooling in the postshock gas can amplify the effect (Sutherland et al. 2003). 

The source of the instability is the use of ID Riemann solvers to compute fluxes in a multidimensional flow. When a planar shock 
is grid aligned, there is too little dissipation added to the fluxes in directions perpendicular to the shock front. Thus, small amplitude 
perturbations in the transverse direction grow, rather than being damped. The solution is to identify grid-aligned shocks and add 
extra dissipation to the transverse fluxes (e.g. Sutherland et al. 2003). In Athena, we use one form of the "H-correction" technique 
described in Sanders et al. (1998) to identify shocks and to add the appropriate dissipation. 

The H-correction is most easily described when used in combination with the Roe fluxes. Consider the calculation of the flux at 
the interface located at (z- 1 /2,j) in 2D. When the H-correction is used the absolute value of the eigenvalues |A"| in the Roe flux 
formula (equation|57|i are replaced with |A a |, where for each component a 

|A Q | = max(|A Q U-_ 1/2J ). (CI) 

Note the max is taken over each |A Q | independently in a pairwise fashion with fji-i/2j, rather than over all a eigenvalues at once. 
Here, ?7/-i/2j comes from a 2D average using a five-point stencil in the shape of the letter 'H', that is 

Vi-l/2,j = ™*(Vi-lJ+l/2,Vi-lJ-l/2,Vi-l/2J,ViJ+l/2,rii.j-l/2) (C2) 

where f?;-i/2j = 5 1 («/ j + C/j j) — (k,-i j — Cf t i-\j)\, Ujj is the component of the velocity normal to the interface, and C/jj is the fast 
magnetosonic speed (for MHD) in the direction normal to the interface. This correction is only added to the final multidimensional 
fluxes (computed in step 6 in 2D, and step 7 in 3D). It only becomes important in shocks, and for grid-aligned shocks it results in the 
dissipation in the transverse directions being comparable to that added in the direction of shock propagation. In 3D the H-correction 
generalizes to a 9-point average (one 'H' in each orthogonal plane). We find the HLL-type fluxes are less susceptible to the carbuncle 
instability, but are still affected by it in some circumstances. The H-correction can be added to HLL-type solvers by making the 
appropriate modification to the wavespeeds b + and b~ defined in equations [53] and [54] 

Use of the H-correction is only required for flows with strong, grid-aligned shocks (for most applications with Athena it is not 
needed). The results of the Noh strong shock test described in jj8.5| show the H-correction is extremely effective at eliminating the 
carbuncle instability. In fact a variety of forms for the correction were proposed by Sanders et al. (see their equation 9). Tests using 
planar shocks in 2D subject to the carbuncle instability showed little difference between the formulations suggested by Sanders et 
al., thus we have chosen to adopt only the version described above. 



